patch bump

move conformance check to avoid doing two checks
simplified fallback in handling of op_sum specialisations
2025-08-08 13:56:43 +10:00 · 2025-08-05 23:51:07 +10:00 · 2025-08-05 23:43:29 +10:00 · 2025-08-05 12:39:40 +10:00 · 2025-07-21 11:51:08 +10:00 · 2025-07-21 11:40:40 +10:00
108 changed files with 3373 additions and 1823 deletions
@@ -16,23 +16,24 @@
 # ------------------------------------------------------------------------


-## NOTE: This installation script by default generates a runtime wrapper library,
+## NOTE: This installation script generates a runtime wrapper library by default,
 ## NOTE: which links with OpenBLAS, LAPACK, etc.
 ## NOTE: 
 ## NOTE: If you prefer to link your programs directly with OpenBLAS, LAPACK, etc
 ## NOTE: there are two options:
+## NOTE: 
 ## NOTE: (1) use the HEADER_ONLY option when running this installation script,
 ## NOTE: or
-## NOTE: (2) do not use this installation script, and instead compile your programs
-## NOTE: using the include folder directly obtained from the armadillo archive.
+## NOTE: (2) do not use this installation script; instead compile your programs
+## NOTE:     using the include folder directly obtained from the armadillo archive.
 ## NOTE: 
 ## NOTE: To forcefully disable use of the runtime wrapper library,
 ## NOTE: compile your programs with ARMA_DONT_USE_WRAPPER enabled.
-## NOTE: For example:
+## NOTE: Example:
 ## NOTE: g++ prog.cpp -o prog -O2 -I /home/user/armadillo-14.0.0/include -DARMA_DONT_USE_WRAPPER -lopenblas -llapack
 ## NOTE: 
-## NOTE: The above command line assumes that you have unpacked the armadillo archive into /home/user/  
-## NOTE: This needs to be adjusted for newer versions of Armadillo (ie. change "14.0.0"),
+## NOTE: The above assumes that the armadillo archive is unpacked into /home/user/  
+## NOTE: Adjust for newer versions of Armadillo (change "14.0.0"),
 ## NOTE: and/or if the armadillo archive has been unpacked into a different folder.
 ## NOTE: 
 ## NOTE: Replace -lopenblas with -lblas if you don't have OpenBLAS.
@@ -40,7 +41,7 @@
 ## NOTE: More details: https://arma.sourceforge.net/faq.html


-cmake_minimum_required(VERSION 3.5 FATAL_ERROR)
+cmake_minimum_required(VERSION 3.5...3.10 FATAL_ERROR)

 if(POLICY CMP0025)
  # enable differentiation between vanilla Clang and AppleClang
@@ -73,7 +74,7 @@ set(ARMA_USE_LAPACK  false)
 set(ARMA_USE_BLAS    false)
 set(ARMA_USE_ATLAS   false)
 set(ARMA_USE_ARPACK  false)
-set(ARMA_USE_SUPERLU false)  # Caveat: only SuperLU version 5.x can be used!
+set(ARMA_USE_SUPERLU false)


 # extract version from sources
@@ -469,7 +470,7 @@ if(ARPACK_FOUND)
  set(ARMA_LIBS ${ARMA_LIBS} ${ARPACK_LIBRARY})
 endif()

-include(ARMA_FindSuperLU5)
+include(ARMA_FindSuperLU)
 message(STATUS "SuperLU_FOUND = ${SuperLU_FOUND}")

 if(SuperLU_FOUND)
@@ -129,8 +129,8 @@ On Windows systems, the MinGW toolset or Visual Studio C++ 2019 (MSVC) can be us

 Caveats on the use of SuperLU:
 - SuperLU must be available as a shared library
- Only the following SuperLU versions are supported: 5.2.x, 5.3.x, 6.0.x
- SuperLU 6.0.x must be compiled with default integer size (32 bits)
+- Only the following SuperLU versions are supported: 5.2.x, 5.3.x, 6.0.x, 7.0.x
+- SuperLU 6.0.x and 7.0.x must be compiled with default integer size (32 bits)

 ---

@@ -369,7 +369,7 @@ Comment out the line containing:

 Use of the C++11 `auto` keyword is not recommended with Armadillo objects and expressions.

-Armadillo has a template meta-programming framework which creates lots of short lived temporaries
+Armadillo has a template meta-programming framework which creates short-lived temporaries
 that are not properly handled by `auto`.

 ---
@@ -0,0 +1,97 @@
+# Try to find suitable version of SuperLU
+# Once done, this will define
+#
+#  SuperLU_FOUND        - system has SuperLU
+#  SuperLU_LIBRARY      - Link this to use SuperLU
+#  SuperLU_INCLUDE_DIR  - directory of SuperLU headers
+
+find_path(SuperLU_INCLUDE_DIR slu_ddefs.h
+  PATHS ${CMAKE_SYSTEM_INCLUDE_PATH} /usr/include /usr/local/include /opt/local/include
+  PATH_SUFFIXES superlu SuperLU ""
+)
+
+find_library(SuperLU_LIBRARY
+  NAMES superlu
+  PATHS ${CMAKE_SYSTEM_LIBRARY_PATH} /usr/lib64 /usr/lib /usr/local/lib64 /usr/local/lib /opt/local/lib64 /opt/local/lib
+)
+
+set(SuperLU_FOUND NO)
+
+if(NOT SuperLU_FIND_QUIETLY)
+  message(STATUS "Looking for suitable version of SuperLU")
+endif()
+
+if(SuperLU_LIBRARY AND SuperLU_INCLUDE_DIR)
+  
+  # determine if this is a suitable version of SuperLU
+  
+  if(EXISTS "${SuperLU_INCLUDE_DIR}/slu_ddefs.h")
+    
+    file(READ "${SuperLU_INCLUDE_DIR}/slu_ddefs.h" _slu_HEADER_CONTENTS)
+    
+    string(REGEX REPLACE ".*version ([0-9]).*" "\\1" SLU_VERSION_MAJOR "${_slu_HEADER_CONTENTS}")
+    
+    if("${SLU_VERSION_MAJOR}" EQUAL "7")
+      set(SuperLU_FOUND YES)
+    endif()
+    
+    if("${SLU_VERSION_MAJOR}" EQUAL "6")
+      set(SuperLU_FOUND YES)
+    endif()
+    
+    if("${SLU_VERSION_MAJOR}" EQUAL "5")
+      set(SuperLU_FOUND YES)
+    endif()
+    
+    if("${SLU_VERSION_MAJOR}" EQUAL "4")
+      
+      # NOTE: SLU_VERSION_MAJOR = 4 is not 100% reliable, as the version string in slu_ddefs.h wasn't updated in SuperLU 5.0, 5.1, 5.2
+      # NOTE: We need to explicitly check the prototype of dgssvx() to determine whether this is version 4 or 5.
+      
+      string(FIND "${_slu_HEADER_CONTENTS}" "dgssvx" _slu_dgssvx_start)
+      
+      if("${_slu_dgssvx_start}" GREATER 0)
+        
+        string(SUBSTRING "${_slu_HEADER_CONTENTS}" ${_slu_dgssvx_start} -1 _slu_dgssvx_tmp)
+        
+        string(FIND "${_slu_dgssvx_tmp}" ";" _slu_dgssvx_len)
+        
+        string(SUBSTRING "${_slu_HEADER_CONTENTS}" ${_slu_dgssvx_start} ${_slu_dgssvx_len} _slu_dgssvx_proto)
+        
+        string(FIND "${_slu_dgssvx_proto}" "GlobalLU_t" _slu_dgssvx_check)
+        
+        if("${_slu_dgssvx_check}" EQUAL -1)
+          # in version 4, dgssvx() doesn't have a GlobalLU_t argument
+          if(NOT SuperLU_FIND_QUIETLY)
+            message(STATUS "Found SuperLU 4, which is insufficient")
+          endif()
+        else()
+          set(SuperLU_FOUND YES)
+        endif()
+      
+      endif()
+    
+    endif()
+  
+  else()
+  
+    message(STATUS "Could not read ${SuperLU_INCLUDE_DIR}/slu_ddefs.h to verify suitable version")
+  
+  endif()
+  
+endif()
+
+if(SuperLU_LIBRARY AND NOT SuperLU_INCLUDE_DIR)
+  if(NOT SuperLU_FIND_QUIETLY)
+    message(STATUS "Incomplete installation of SuperLU: no SuperLU headers found")
+  endif()
+endif()
+  
+if(NOT SuperLU_FOUND) 
+  if(NOT SuperLU_FIND_QUIETLY)
+    message(STATUS "Could not find suitable version of SuperLU")
+  endif()
+else()
+  message(STATUS "Found SuperLU library: ${SuperLU_LIBRARY}")
+endif()
+
@@ -1,90 +0,0 @@
-# Try to find SuperLU 5.x
-# Once done, this will define
-#
-#  SuperLU_FOUND        - system has SuperLU
-#  SuperLU_LIBRARY      - Link this to use SuperLU
-#  SuperLU_INCLUDE_DIR  - directory of SuperLU headers
-
-find_path(SuperLU_INCLUDE_DIR slu_ddefs.h
-  PATHS ${CMAKE_SYSTEM_INCLUDE_PATH} /usr/include /usr/local/include /opt/local/include
-  PATH_SUFFIXES superlu SuperLU ""
-)
-
-find_library(SuperLU_LIBRARY
-  NAMES superlu
-  PATHS ${CMAKE_SYSTEM_LIBRARY_PATH} /usr/lib64 /usr/lib /usr/local/lib64 /usr/local/lib /opt/local/lib64 /opt/local/lib
-)
-
-set(SuperLU_FOUND NO)
-
-if(NOT SuperLU_FIND_QUIETLY)
-  message(STATUS "Looking for SuperLU version 5")
-endif()
-
-if(SuperLU_LIBRARY AND SuperLU_INCLUDE_DIR)
-  
-  # Now determine if the version is version 5.
-  # This voodoo is required as SuperLU developers like to keep things messy.
-  
-  if(EXISTS "${SuperLU_INCLUDE_DIR}/slu_ddefs.h")
-    
-    file(READ "${SuperLU_INCLUDE_DIR}/slu_ddefs.h" _slu_HEADER_CONTENTS)
-    
-    string(REGEX REPLACE ".*version ([0-9]).*" "\\1" SLU_VERSION_MAJOR "${_slu_HEADER_CONTENTS}")
-    
-    if("${SLU_VERSION_MAJOR}" EQUAL "5")
-      set(SuperLU_FOUND YES)
-    else()
-      if("${SLU_VERSION_MAJOR}" EQUAL "4")
-        
-        # NOTE: SLU_VERSION_MAJOR = 4 is not 100% reliable, as the version string in slu_ddefs.h wasn't updated in SuperLU 5.0, 5.1, 5.2
-        # NOTE: We need to explicitly check the prototype of dgssvx() to determine whether this is version 4 or 5.
-        
-        string(FIND "${_slu_HEADER_CONTENTS}" "dgssvx" _slu_dgssvx_start)
-        
-        if("${_slu_dgssvx_start}" GREATER 0)
-          
-          string(SUBSTRING "${_slu_HEADER_CONTENTS}" ${_slu_dgssvx_start} -1 _slu_dgssvx_tmp)
-          
-          string(FIND "${_slu_dgssvx_tmp}" ";" _slu_dgssvx_len)
-          
-          string(SUBSTRING "${_slu_HEADER_CONTENTS}" ${_slu_dgssvx_start} ${_slu_dgssvx_len} _slu_dgssvx_proto)
-          
-          string(FIND "${_slu_dgssvx_proto}" "GlobalLU_t" _slu_dgssvx_check)
-          
-          if("${_slu_dgssvx_check}" EQUAL -1)
-            # in version 4, dgssvx() doesn't have a GlobalLU_t argument
-            if(NOT SuperLU_FIND_QUIETLY)
-              message(STATUS "Found SuperLU, but it doesn't appear to be version 5")
-            endif()
-          else()
-            set(SuperLU_FOUND YES)
-          endif()
-        endif()
-      else()
-        if(NOT SuperLU_FIND_QUIETLY)
-          message(STATUS "Found SuperLU version ${SLU_VERSION_MAJOR}, but version 5 is required")
-        endif()
-      endif()
-    endif()
-  else()
-    if(NOT SuperLU_FIND_QUIETLY)
-      message(STATUS "Found SuperLU and headers, but could not verify version 5")
-    endif()
-  endif()
-endif()
-
-if(SuperLU_LIBRARY AND NOT SuperLU_INCLUDE_DIR)
-  if(NOT SuperLU_FIND_QUIETLY)
-    message(STATUS "Incomplete installation of SuperLU: no SuperLU headers found")
-  endif()
-endif()
-  
-if(NOT SuperLU_FOUND) 
-  if(NOT SuperLU_FIND_QUIETLY)
-    message(STATUS "Could not find SuperLU")
-  endif()
-else()
-  message(STATUS "Found SuperLU: ${SuperLU_LIBRARY}")
-endif()
-
@@ -304,6 +304,7 @@ namespace arma
  #include "armadillo_bits/op_clamp_bones.hpp"
  #include "armadillo_bits/op_expmat_bones.hpp"
  #include "armadillo_bits/op_nonzeros_bones.hpp"
+  #include "armadillo_bits/op_omit_bones.hpp"
  #include "armadillo_bits/op_diff_bones.hpp"
  #include "armadillo_bits/op_norm_bones.hpp"
  #include "armadillo_bits/op_vecnorm_bones.hpp"
@@ -377,6 +378,7 @@ namespace arma
  #include "armadillo_bits/spop_norm_bones.hpp"
  #include "armadillo_bits/spop_shift_bones.hpp"
  #include "armadillo_bits/spop_relational_bones.hpp"
+  #include "armadillo_bits/spop_omit_bones.hpp"
  
  #include "armadillo_bits/spglue_plus_bones.hpp"
  #include "armadillo_bits/spglue_minus_bones.hpp"
@@ -564,6 +566,7 @@ namespace arma
  #include "armadillo_bits/fn_clamp.hpp"
  #include "armadillo_bits/fn_expmat.hpp"
  #include "armadillo_bits/fn_nonzeros.hpp"
+  #include "armadillo_bits/fn_omit.hpp"
  #include "armadillo_bits/fn_interp1.hpp"
  #include "armadillo_bits/fn_interp2.hpp"
  #include "armadillo_bits/fn_qz.hpp"
@@ -590,6 +593,7 @@ namespace arma
  #include "armadillo_bits/fn_powmat.hpp"
  #include "armadillo_bits/fn_powext.hpp"
  #include "armadillo_bits/fn_diags_spdiags.hpp"
+  #include "armadillo_bits/fn_balance.hpp"
  
  #include "armadillo_bits/fn_speye.hpp"
  #include "armadillo_bits/fn_spones.hpp"
@@ -758,6 +762,7 @@ namespace arma
  #include "armadillo_bits/op_clamp_meat.hpp"
  #include "armadillo_bits/op_expmat_meat.hpp"
  #include "armadillo_bits/op_nonzeros_meat.hpp"
+  #include "armadillo_bits/op_omit_meat.hpp"
  #include "armadillo_bits/op_diff_meat.hpp"
  #include "armadillo_bits/op_norm_meat.hpp"
  #include "armadillo_bits/op_vecnorm_meat.hpp"
@@ -831,6 +836,7 @@ namespace arma
  #include "armadillo_bits/spop_norm_meat.hpp"
  #include "armadillo_bits/spop_shift_meat.hpp"
  #include "armadillo_bits/spop_relational_meat.hpp"
+  #include "armadillo_bits/spop_omit_meat.hpp"
  
  #include "armadillo_bits/spglue_plus_meat.hpp"
  #include "armadillo_bits/spglue_minus_meat.hpp"
@@ -328,7 +328,7 @@ BaseCube<elem_type,derived>::is_finite() const
    for(uword c=0; c<n_c; ++c)
    for(uword r=0; r<n_r; ++r)
      {
-      if( arma_isfinite(P.at(r,c,s)) == false )  { return false; }
+      if( arma_isnonfinite(P.at(r,c,s)) )  { return false; }
      }
    }
  
@@ -436,7 +436,7 @@ BaseCube<elem_type,derived>::has_nonfinite() const
    for(uword c=0; c<n_c; ++c)
    for(uword r=0; r<n_r; ++r)
      {
-      if(arma_isfinite(P.at(r,c,s)) == false)  { return true; }
+      if(arma_isnonfinite(P.at(r,c,s)))  { return true; }
      }
    }
  
@@ -685,7 +685,7 @@ Base<elem_type,derived>::is_finite() const
      
      for(uword i=0; i<n_elem; ++i)
        {
-        if(arma_isfinite(Pea[i]) == false)  { return false; }
+        if(arma_isnonfinite(Pea[i]))  { return false; }
        }
      }
    else
@@ -696,7 +696,7 @@ Base<elem_type,derived>::is_finite() const
      for(uword col=0; col<n_cols; ++col)
      for(uword row=0; row<n_rows; ++row)
        {
-        if(arma_isfinite(P.at(row,col)) == false)  { return false; }
+        if(arma_isnonfinite(P.at(row,col)))  { return false; }
        }
      }
    }
@@ -829,7 +829,7 @@ Base<elem_type,derived>::has_nonfinite() const
      
      for(uword i=0; i<n_elem; ++i)
        {
-        if(arma_isfinite(Pea[i]) == false)  { return true; }
+        if(arma_isnonfinite(Pea[i]))  { return true; }
        }
      }
    else
@@ -840,7 +840,7 @@ Base<elem_type,derived>::has_nonfinite() const
      for(uword col=0; col<n_cols; ++col)
      for(uword row=0; row<n_rows; ++row)
        {
-        if(arma_isfinite(P.at(row,col)) == false)  { return true; }
+        if(arma_isnonfinite(P.at(row,col)))  { return true; }
        }
      }
    }
@@ -784,6 +784,10 @@ Cube<eT>::operator=(const Cube<eT>& x)
    
    arrayops::copy( memptr(), x.mem, n_elem );
    }
+  else
+    {
+    arma_debug_print("Cube::operator=(): copy omitted");
+    }
  
  return *this;
  }
@@ -2868,7 +2872,7 @@ Cube<eT>::Cube(const eOpCube<T1, eop_type>& X)
  
  init_cold();
  
-  if(is_same_type<eop_type, eop_pow>::value)
+  if(arma_config::optimise_powexpr && is_same_type<eop_type, eop_pow>::value)
    {
    constexpr bool eT_non_int = is_non_integral<eT>::value;
    
@@ -2898,7 +2902,7 @@ Cube<eT>::operator=(const eOpCube<T1, eop_type>& X)
  
  init_warm(X.get_n_rows(), X.get_n_cols(), X.get_n_slices());
  
-  if(is_same_type<eop_type, eop_pow>::value)
+  if(arma_config::optimise_powexpr && is_same_type<eop_type, eop_pow>::value)
    {
    constexpr bool eT_non_int = is_non_integral<eT>::value;
    
@@ -2928,7 +2932,7 @@ Cube<eT>::operator+=(const eOpCube<T1, eop_type>& X)
  
  if(bad_alias)  { const Cube<eT> tmp(X); return (*this).operator+=(tmp); }
  
-  if(is_same_type<eop_type, eop_pow>::value)
+  if(arma_config::optimise_powexpr && is_same_type<eop_type, eop_pow>::value)
    {
    constexpr bool eT_non_int = is_non_integral<eT>::value;
    
@@ -2958,7 +2962,7 @@ Cube<eT>::operator-=(const eOpCube<T1, eop_type>& X)
  
  if(bad_alias)  { const Cube<eT> tmp(X); return (*this).operator-=(tmp); }
  
-  if(is_same_type<eop_type, eop_pow>::value)
+  if(arma_config::optimise_powexpr && is_same_type<eop_type, eop_pow>::value)
    {
    constexpr bool eT_non_int = is_non_integral<eT>::value;
    
@@ -2988,7 +2992,7 @@ Cube<eT>::operator%=(const eOpCube<T1, eop_type>& X)
  
  if(bad_alias)  { const Cube<eT> tmp(X); return (*this).operator%=(tmp); }
  
-  if(is_same_type<eop_type, eop_pow>::value)
+  if(arma_config::optimise_powexpr && is_same_type<eop_type, eop_pow>::value)
    {
    constexpr bool eT_non_int = is_non_integral<eT>::value;
    
@@ -3018,7 +3022,7 @@ Cube<eT>::operator/=(const eOpCube<T1, eop_type>& X)
  
  if(bad_alias)  { const Cube<eT> tmp(X); return (*this).operator/=(tmp); }
  
-  if(is_same_type<eop_type, eop_pow>::value)
+  if(arma_config::optimise_powexpr && is_same_type<eop_type, eop_pow>::value)
    {
    constexpr bool eT_non_int = is_non_integral<eT>::value;
    
@@ -106,7 +106,12 @@ MapMat<eT>::operator=(const MapMat<eT>& x)
  {
  arma_debug_sigprint();
  
-  if(this == &x)  { return; }
+  if(this == &x)
+    {
+    arma_debug_print("MapMat::operator=(): copy omitted");
+    
+    return;
+    }
  
  access::rw(n_rows) = x.n_rows;
  access::rw(n_cols) = x.n_cols;
@@ -970,18 +975,17 @@ MapMat_val<eT>::operator*=(const eT in_val)
  
  if(it != it_end)
    {
-    if(in_val != eT(0))
-      {
-      eT& val = (*it).second;
-      
-      val *= in_val;
-      
-      if(val == eT(0))  { map_ref.erase(it); }
-      }
-    else
-      {
-      map_ref.erase(it);
-      }
+    eT& val = (*it).second;
+    
+    val *= in_val;
+    
+    if(val == eT(0))  { map_ref.erase(it); }
+    }
+  else
+    {
+    const eT val = eT(0) * in_val;  // in case in_val is inf or nan
+    
+    if(val != eT(0))  { parent.set_val(index, val); }
    }
  }

@@ -1009,9 +1013,7 @@ MapMat_val<eT>::operator/=(const eT in_val)
    }
  else
    {
-    // silly operation, but included for completeness
-    
-    const eT val = eT(0) / in_val;
+    const eT val = eT(0) / in_val;  // in case in_val is zero or nan
    
    if(val != eT(0))  { parent.set_val(index, val); }
    }
@@ -1481,18 +1483,11 @@ SpMat_MapMat_val<eT>::mul(const eT in_val)
    
    if(it != it_end)
      {
-      if(in_val != eT(0))
-        {
-        eT& val = (*it).second;
-        
-        val *= in_val;
-        
-        if(val == eT(0))  { map_ref.erase(it); }
-        }
-      else
-        {
-        map_ref.erase(it);
-        }
+      eT& val = (*it).second;
+      
+      val *= in_val;
+      
+      if(val == eT(0))  { map_ref.erase(it); }
      
      s_parent.sync_state = 1;
      
@@ -1500,19 +1495,15 @@ SpMat_MapMat_val<eT>::mul(const eT in_val)
      }
    else
      {
-      // element not found, ie. it's zero; zero multiplied by anything is zero, except for nan and inf
-      if(arma_isfinite(in_val) == false)
+      const eT result = eT(0) * in_val;  // in case in_val is inf or nan
+      
+      if(result != eT(0))
        {
-        const eT result = eT(0) * in_val;
+        m_parent.set_val(index, result);
        
-        if(result != eT(0))  // paranoia, in case compiling with -ffast-math
-          {
-          m_parent.set_val(index, result);
-          
-          s_parent.sync_state = 1;
-          
-          access::rw(s_parent.n_nonzero) = m_parent.get_n_nonzero();
-          }
+        s_parent.sync_state = 1;
+        
+        access::rw(s_parent.n_nonzero) = m_parent.get_n_nonzero();
        }
      }
    }
@@ -1554,19 +1545,15 @@ SpMat_MapMat_val<eT>::div(const eT in_val)
      }
    else
      {
-      // element not found, ie. it's zero; zero divided by anything is zero, except for zero and nan
-      if( (in_val == eT(0)) || (arma_isnan(in_val)) )
+      const eT result = eT(0) / in_val;  // in case in_val is zero or nan
+      
+      if(result != eT(0))
        {
-        const eT result = eT(0) / in_val;
+        m_parent.set_val(index, result);
        
-        if(result != eT(0))  // paranoia, in case compiling with -ffast-math
-          {
-          m_parent.set_val(index, result);
-          
-          s_parent.sync_state = 1;
-          
-          access::rw(s_parent.n_nonzero) = m_parent.get_n_nonzero();
-          }
+        s_parent.sync_state = 1;
+        
+        access::rw(s_parent.n_nonzero) = m_parent.get_n_nonzero();
        }
      }
    }
@@ -915,6 +915,10 @@ Mat<eT>::operator=(const Mat<eT>& in_mat)
    
    arrayops::copy( memptr(), in_mat.mem, in_mat.n_elem );
    }
+  else
+    {
+    arma_debug_print("Mat::operator=(): copy omitted");
+    }
  
  return *this;
  }
@@ -5197,7 +5201,7 @@ Mat<eT>::Mat(const eOp<T1, eop_type>& X)
  
  init_cold();
  
-  if(is_same_type<eop_type, eop_pow>::value)
+  if(arma_config::optimise_powexpr && is_same_type<eop_type, eop_pow>::value)
    {
    constexpr bool eT_non_int = is_non_integral<eT>::value;
    
@@ -5227,7 +5231,7 @@ Mat<eT>::operator=(const eOp<T1, eop_type>& X)
  
  init_warm(X.get_n_rows(), X.get_n_cols());
  
-  if(is_same_type<eop_type, eop_pow>::value)
+  if(arma_config::optimise_powexpr && is_same_type<eop_type, eop_pow>::value)
    {
    constexpr bool eT_non_int = is_non_integral<eT>::value;
    
@@ -5256,7 +5260,7 @@ Mat<eT>::operator+=(const eOp<T1, eop_type>& X)
  
  if(bad_alias)  { const Mat<eT> tmp(X); return (*this).operator+=(tmp); }
  
-  if(is_same_type<eop_type, eop_pow>::value)
+  if(arma_config::optimise_powexpr && is_same_type<eop_type, eop_pow>::value)
    {
    constexpr bool eT_non_int = is_non_integral<eT>::value;
    
@@ -5285,7 +5289,7 @@ Mat<eT>::operator-=(const eOp<T1, eop_type>& X)
  
  if(bad_alias)  { const Mat<eT> tmp(X); return (*this).operator-=(tmp); }
  
-  if(is_same_type<eop_type, eop_pow>::value)
+  if(arma_config::optimise_powexpr && is_same_type<eop_type, eop_pow>::value)
    {
    constexpr bool eT_non_int = is_non_integral<eT>::value;
    
@@ -5331,7 +5335,7 @@ Mat<eT>::operator%=(const eOp<T1, eop_type>& X)
  
  if(bad_alias)  { const Mat<eT> tmp(X); return (*this).operator%=(tmp); }
  
-  if(is_same_type<eop_type, eop_pow>::value)
+  if(arma_config::optimise_powexpr && is_same_type<eop_type, eop_pow>::value)
    {
    constexpr bool eT_non_int = is_non_integral<eT>::value;
    
@@ -5360,7 +5364,7 @@ Mat<eT>::operator/=(const eOp<T1, eop_type>& X)
  
  if(bad_alias)  { const Mat<eT> tmp(X); return (*this).operator/=(tmp); }
  
-  if(is_same_type<eop_type, eop_pow>::value)
+  if(arma_config::optimise_powexpr && is_same_type<eop_type, eop_pow>::value)
    {
    constexpr bool eT_non_int = is_non_integral<eT>::value;
    
@@ -748,7 +748,7 @@ SpBase<elem_type,derived>::is_finite() const
    
    while(it != it_end)
      {
-      if(arma_isfinite(*it) == false)  { return false; }
+      if(arma_isnonfinite(*it))  { return false; }
      ++it;
      }
    }
@@ -850,7 +850,7 @@ SpBase<elem_type,derived>::has_nonfinite() const
    
    while(it != it_end)
      {
-      if(arma_isfinite(*it) == false)  { return true; }
+      if(arma_isnonfinite(*it))  { return true; }
      ++it;
      }
    }
@@ -612,33 +612,26 @@ SpMat<eT>::operator*=(const eT val)
  {
  arma_debug_sigprint();
  
-  if(val != eT(0))
+  sync_csc();
+  invalidate_cache();
+  
+  const uword n_nz = n_nonzero;
+  
+  eT* vals = access::rwp(values);
+  
+  bool has_zero = false;
+  
+  for(uword i=0; i<n_nz; ++i)
    {
-    sync_csc();
-    invalidate_cache();
+    eT& vals_i = vals[i];
    
-    const uword n_nz = n_nonzero;
+    vals_i *= val;
    
-    eT* vals = access::rwp(values);
-    
-    bool has_zero = false;
-    
-    for(uword i=0; i<n_nz; ++i)
-      {
-      eT& vals_i = vals[i];
-      
-      vals_i *= val;
-      
-      if(vals_i == eT(0))  { has_zero = true; }
-      }
-    
-    if(has_zero)  { remove_zeros(); }
-    }
-  else
-    {
-    (*this).zeros();
+    if(vals_i == eT(0))  { has_zero = true; }
    }
  
+  if(has_zero)  { remove_zeros(); }
+  
  return *this;
  }

@@ -5131,7 +5124,12 @@ SpMat<eT>::init(const SpMat<eT>& x)
  {
  arma_debug_sigprint();
  
-  if(this == &x)  { return; }
+  if(this == &x)
+    {
+    arma_debug_print("SpMat::init(): copy omitted");
+    
+    return;
+    }
  
  bool init_done = false;
  
@@ -253,8 +253,6 @@ SpSubview_col_list<eT,T1>::operator*= (const eT val)
  {
  arma_debug_sigprint();
  
-  if(val == eT(0))  { (*this).zeros(); return; }
-  
  SpMat<eT>& m_local = const_cast< SpMat<eT>& >(m);
  
  const umat&  ci        = U_ci.M;
@@ -157,8 +157,6 @@ SpSubview<eT>::operator*=(const eT val)
  {
  arma_debug_sigprint();
  
-  if(val == eT(0))  { (*this).zeros(); return *this; }
-  
  if((n_elem == 0) || (n_nonzero == 0))  { return *this; }
  
  m.sync_csc();
@@ -148,24 +148,17 @@ inline
 SpValProxy<T1>&
 SpValProxy<T1>::operator*=(const eT rhs)
  {
-  if(rhs != eT(0))
+  if(val_ptr)
    {
-    if(val_ptr)
-      {
-      // The value already exists and merely needs to be updated.
-      *val_ptr *= rhs;
-      parent.invalidate_cache();
-      check_zero();
-      }
+    *val_ptr *= rhs;
+    parent.invalidate_cache();
+    check_zero();
    }
  else
    {
-    if(val_ptr)
-      {
-      // Since we are multiplying by zero, the value can be deleted.
-      parent.delete_element(row, col);
-      val_ptr = nullptr;
-      }
+    const eT val = eT(0) * rhs;  // in case rhs is inf or nan
+    
+    if(val != eT(0))  { val_ptr = &parent.insert_element(row, col, val); }
    }
  
  return *this;
@@ -178,37 +171,17 @@ inline
 SpValProxy<T1>&
 SpValProxy<T1>::operator/=(const eT rhs)
  {
-  if(rhs != eT(0)) // I hope this is true!
+  if(val_ptr)
    {
-    if(val_ptr)
-      {
-      *val_ptr /= rhs;
-      parent.invalidate_cache();
-      check_zero();
-      }
+    *val_ptr /= rhs;
+    parent.invalidate_cache();
+    check_zero();
    }
  else
    {
-    if(val_ptr)
-      {
-      *val_ptr /= rhs; // That is where it gets ugly.
-      // Now check if it's 0.
-      if(*val_ptr == eT(0))
-        {
-        parent.delete_element(row, col);
-        val_ptr = nullptr;
-        }
-      }
-    else
-      {
-      eT val = eT(0) / rhs; // This may vary depending on type and implementation.
-      
-      if(val != eT(0))
-        {
-        // Ok, now we have to insert it.
-        val_ptr = &parent.insert_element(row, col, val);
-        }
-      }
+    const eT val = eT(0) / rhs;  // in case rhs is zero or nan
+    
+    if(val != eT(0))  { val_ptr = &parent.insert_element(row, col, val); }
    }
  
  return *this;
@@ -65,6 +65,48 @@ arma_isfinite(const std::complex<T>& x)
  }


+//
+
+
+template<typename eT>
+inline
+bool
+arma_isnonfinite(eT)
+  {
+  return false;
+  }
+
+
+
+template<>
+inline
+bool
+arma_isnonfinite(float x)
+  {
+  return (std::isfinite(x) == false);
+  }
+
+
+
+template<>
+inline
+bool
+arma_isnonfinite(double x)
+  {
+  return (std::isfinite(x) == false);
+  }
+
+
+
+template<typename T>
+inline
+bool
+arma_isnonfinite(const std::complex<T>& x)
+  {
+  return ( (std::isfinite(x.real()) == false) || (std::isfinite(x.imag()) == false) );
+  }
+
+

 //
 // wrappers for isinf
@@ -65,6 +65,13 @@ struct arma_config
  #endif
  
  
+  #if defined(ARMA_OPTIMISE_POWEXPR)
+    static constexpr bool optimise_powexpr = true;
+  #else
+    static constexpr bool optimise_powexpr = false;
+  #endif
+  
+  
  #if defined(ARMA_CHECK_CONFORMANCE)
    static constexpr bool check_conform = true;
  #else
@@ -93,6 +93,7 @@ class op_vectorise_row;
 class op_vectorise_col;
 class op_symmatu;
 class op_symmatl;
+class op_omit;

 class op_row_as_mat;
 class op_col_as_mat;
@@ -361,6 +362,7 @@ struct arma_zeros_indicator   : public arma_initmode_indicator<true > {};
 struct arma_nozeros_indicator : public arma_initmode_indicator<false> {};


+
 //! \addtogroup injector
 //! @{

@@ -70,7 +70,7 @@ arma_ostream::modify_stream(std::ostream& o, const eT* data, const uword n_elem)
    {
    const eT val = data[i];
    
-    if(arma_isfinite(val) == false)  { continue; }
+    if(arma_isnonfinite(val))  { continue; }
    
    if(
      ( cond_rel< (sizeof(eT) > 4) && (is_same_type<uword,eT>::yes || is_same_type<sword,eT>::yes) >::geq(val, eT(+10000000000)) )
@@ -208,7 +208,7 @@ arma_ostream::modify_stream(std::ostream& o, typename SpMat<eT>::const_iterator
    {
    const eT val = (*it);
    
-    if(arma_isfinite(val) == false)  { continue; }
+    if(arma_isnonfinite(val))  { continue; }
    
    if(
      val >= eT(+100) ||
@@ -22,9 +22,9 @@


 #define ARMA_VERSION_MAJOR 14
-#define ARMA_VERSION_MINOR 4
-#define ARMA_VERSION_PATCH 1
-#define ARMA_VERSION_NAME  "Filtered Espresso"
+#define ARMA_VERSION_MINOR 6
+#define ARMA_VERSION_PATCH 2
+#define ARMA_VERSION_NAME  "Caffe Mocha"



@@ -1050,13 +1050,13 @@ arrayops::is_finite(const eT* src, const uword n_elem)
    const eT val_i = (*src);  src++;
    const eT val_j = (*src);  src++;
    
-    if(arma_isfinite(val_i) == false)  { return false; }
-    if(arma_isfinite(val_j) == false)  { return false; }
+    if(arma_isnonfinite(val_i))  { return false; }
+    if(arma_isnonfinite(val_j))  { return false; }
    }
  
  if((j-1) < n_elem)
    {
-    if(arma_isfinite(*src) == false)  { return false; }
+    if(arma_isnonfinite(*src))  { return false; }
    }
  
  return true;
@@ -391,7 +391,7 @@ class auxlib
  // solve the Sylvester equation AX + XB = C
  
  template<typename eT>
-  inline static bool syl(Mat<eT>& X, const Mat<eT>& A, const Mat<eT>& B, const Mat<eT>& C);
+  inline static bool sylvester(Mat<eT>& X, const Mat<eT>& A, const Mat<eT>& B, const Mat<eT>& C);
  
  
  //
@@ -404,6 +404,13 @@ class auxlib
  inline static bool qz(Mat< std::complex<T> >& A, Mat< std::complex<T> >& B, Mat< std::complex<T> >& vsl, Mat< std::complex<T> >& vsr, const Base< std::complex<T>, T1 >& X_expr, const Base< std::complex<T>, T2 >& Y_expr, const char mode);
  
  
+  //
+  // matrix balance
+  
+  template<typename eT>
+  inline static bool balance(Col<typename get_pod_type<eT>::result>& S, Col<uword>& P, Mat<eT>& A, const bool calc_SP, const bool do_scal, const bool do_perm);
+  
+  
  // 
  // rcond
  
@@ -6486,15 +6486,15 @@ auxlib::schur(Mat< std::complex<T> >& U, Mat< std::complex<T> >& S, const bool c
 template<typename eT>
 inline
 bool
-auxlib::syl(Mat<eT>& X, const Mat<eT>& A, const Mat<eT>& B, const Mat<eT>& C)
+auxlib::sylvester(Mat<eT>& X, const Mat<eT>& A, const Mat<eT>& B, const Mat<eT>& C)
  {
  arma_debug_sigprint();
  
  #if defined(ARMA_USE_LAPACK)
    {
-    arma_conform_check( (A.is_square() == false) || (B.is_square() == false), "syl(): given matrices must be square sized" );
+    arma_conform_check( (A.is_square() == false) || (B.is_square() == false), "sylvester(): given matrices must be square sized" );
      
-    arma_conform_check( (C.n_rows != A.n_rows) || (C.n_cols != B.n_cols), "syl(): matrices are not conformant" );
+    arma_conform_check( (C.n_rows != A.n_rows) || (C.n_cols != B.n_cols), "sylvester(): matrices are not conformant" );
    
    if(A.is_empty() || B.is_empty() || C.is_empty())  { X.reset(); return true; }
    
@@ -6534,7 +6534,7 @@ auxlib::syl(Mat<eT>& X, const Mat<eT>& A, const Mat<eT>& B, const Mat<eT>& C)
    arma_ignore(A);
    arma_ignore(B);
    arma_ignore(C);
-    arma_stop_logic_error("syl(): use of LAPACK must be enabled");
+    arma_stop_logic_error("sylvester(): use of LAPACK must be enabled");
    return false;
    }
  #endif
@@ -6714,6 +6714,76 @@ auxlib::qz(Mat< std::complex<T> >& A, Mat< std::complex<T> >& B, Mat< std::compl



+template<typename eT>
+inline
+bool
+auxlib::balance(Col<typename get_pod_type<eT>::result>& S, Col<uword>& P, Mat<eT>& A, const bool calc_SP, const bool do_scal, const bool do_perm)
+  {
+  arma_debug_sigprint();
+  
+  #if defined(ARMA_USE_LAPACK)
+    {
+    typedef typename get_pod_type<eT>::result T;
+    
+    // assuming given matrix is square-sized
+    
+    if(A.n_elem == 0)  { S.reset(); P.reset(); return true; }
+    
+    const char job = (do_scal && do_perm) ? 'B' : ((do_scal) ? 'S' : ((do_perm) ? 'P' : 'N'));
+    
+    blas_int n    = blas_int(A.n_rows);
+    blas_int lda  = blas_int(A.n_rows);
+    blas_int ilo  = blas_int(0);
+    blas_int ihi  = blas_int(0);
+    blas_int info = blas_int(0);
+    
+    podarray<T> scale(A.n_rows);  scale.zeros();
+    
+    arma_debug_print("lapack::gebal()");
+    lapack::gebal(&job, &n, A.memptr(), &lda, &ilo, &ihi, scale.memptr(), &info);
+    
+    if(info != blas_int(0))  { return false; }
+    
+    if(calc_SP == false)  { return true; }
+    
+    const uword N = A.n_rows;
+    
+    // sanity check
+    if( (ilo < 1) || (uword(ihi) > N) )  { arma_debug_print("ilo and/or ihi out of bounds"); return false; }
+    
+    S.zeros(N);
+    P.zeros(N);
+    
+        T* S_mem = S.memptr();
+    uword* P_mem = P.memptr();
+    
+    const T* scale_mem = scale.memptr();
+    
+    for(uword i = 0;            i < uword(ilo)-1; ++i)  { S_mem[i] = T(1);         }
+    for(uword i = uword(ilo)-1; i < uword(ihi);   ++i)  { S_mem[i] = scale_mem[i]; }
+    for(uword i = uword(ihi);   i < N;            ++i)  { S_mem[i] = T(1);         }
+    
+    for(uword i=0; i < N; ++i)  { P_mem[i] = i; }
+    
+    for(uword i=N-1; i >= uword(ihi)  ; --i)  { const uword j = uword(scale_mem[i]) - 1; std::swap(P_mem[i], P_mem[j]); }
+    for(uword i=0;   i <  uword(ilo)-1; ++i)  { const uword j = uword(scale_mem[i]) - 1; std::swap(P_mem[i], P_mem[j]); }
+    
+    return true;
+    }
+  #else
+    {
+    arma_ignore(S);
+    arma_ignore(P);
+    arma_ignore(A);
+    arma_ignore(do_scal);
+    arma_ignore(do_perm);
+    return false;
+    }
+  #endif
+  }
+
+
+
 template<typename eT>
 inline
 eT
@@ -169,7 +169,7 @@
  // https://gcc.gnu.org/bugzilla/show_bug.cgi?id=57580
  
  #if (ARMA_GCC_VERSION < 80100)
-    #pragma message("INFO: support for GCC versions older than 8.1 is deprecated"
+    #pragma message("INFO: support for GCC versions older than 8.1 is deprecated")
  #endif
  
  #define ARMA_GOOD_COMPILER
@@ -477,9 +477,6 @@
 #undef major


-// WARNING: option 'ARMA_IGNORE_DEPRECATED_MARKER' will be removed;
-// WARNING: disabling deprecation messages is counter-productive
-
 #if defined(ARMA_IGNORE_DEPRECATED_MARKER)
  #undef  arma_deprecated
  #define arma_deprecated
@@ -166,6 +166,11 @@
  //// of inv() and inv_sympd() within compound expressions
 #endif

+#if !defined(ARMA_OPTIMISE_POWEXPR)
+  #define ARMA_OPTIMISE_POWEXPR
+  //// Comment out the above line to disable optimised handling of pow()
+#endif
+
 #if !defined(ARMA_CHECK_CONFORMANCE)
  #define ARMA_CHECK_CONFORMANCE
  //// Comment out the above line to disable conformance checks for bounds and size.
@@ -211,7 +216,7 @@


 #if defined(ARMA_DEFAULT_OSTREAM)
-  #pragma message ("WARNING: support for ARMA_DEFAULT_OSTREAM is deprecated and will be removed;")
+  #pragma message ("WARNING: option ARMA_DEFAULT_OSTREAM is deprecated and will be removed;")
  #pragma message ("WARNING: use ARMA_COUT_STREAM and ARMA_CERR_STREAM instead")
 #endif

@@ -290,7 +295,7 @@

 // for compatibility with earlier versions of Armadillo
 #if defined(ARMA_DONT_USE_CXX11_MUTEX)
-  #pragma message ("WARNING: support for ARMA_DONT_USE_CXX11_MUTEX is deprecated and will be removed;")
+  #pragma message ("WARNING: option ARMA_DONT_USE_CXX11_MUTEX is deprecated and will be removed;")
  #pragma message ("WARNING: use ARMA_DONT_USE_STD_MUTEX instead")
  #undef ARMA_USE_STD_MUTEX
 #endif
@@ -303,6 +308,18 @@
  #undef ARMA_64BIT_WORD
 #endif

+#if (defined(ARMA_BLAS_LONG_LONG) && defined(ARMA_USE_WRAPPER))
+  #pragma message ("WARNING: use of ARMA_BLAS_LONG_LONG in conjunction with ARMA_USE_WRAPPER is not supported")
+#endif
+
+#if (defined(ARMA_BLAS_64BIT_INT) && defined(ARMA_USE_WRAPPER))
+  #pragma message ("WARNING: use of ARMA_BLAS_64BIT_INT in conjunction with ARMA_USE_WRAPPER is not supported")
+#endif
+
+#if (defined(ARMA_SUPERLU_64BIT_INT) && defined(ARMA_USE_WRAPPER))
+  #pragma message ("WARNING: use of ARMA_SUPERLU_64BIT_INT in conjunction with ARMA_USE_WRAPPER is not supported")
+#endif
+
 // for compatibility with earlier versions of Armadillo
 #if defined(ARMA_BLAS_LONG) || defined(ARMA_BLAS_LONG_LONG)
  #undef  ARMA_BLAS_64BIT_INT
@@ -321,6 +338,10 @@
  #undef ARMA_OPTIMISE_INVEXPR
 #endif

+#if defined(ARMA_DONT_OPTIMISE_POWEXPR)
+  #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")
@@ -366,12 +387,12 @@
  #undef ARMA_CRIPPLED_LAPACK
 #endif

-// WARNING: option ARMA_IGNORE_DEPRECATED_MARKER will be removed
-// WARNING: option ARMA_CRIPPLED_LAPACK          will be removed
+// WARNING: option ARMA_IGNORE_DEPRECATED_MARKER is deprecated and will be removed
+// WARNING: option ARMA_CRIPPLED_LAPACK          is deprecated and will be removed

 #if defined(ARMA_CRIPPLED_LAPACK)
  #if (!defined(ARMA_IGNORE_DEPRECATED_MARKER))
-    #pragma message ("option ARMA_CRIPPLED_LAPACK is deprecated and will be removed")
+    #pragma message ("WARNING: option ARMA_CRIPPLED_LAPACK is deprecated and will be removed")
  #endif
 #endif

@@ -166,6 +166,11 @@
  //// of inv() and inv_sympd() within compound expressions
 #endif

+#if !defined(ARMA_OPTIMISE_POWEXPR)
+  #define ARMA_OPTIMISE_POWEXPR
+  //// Comment out the above line to disable optimised handling of pow()
+#endif
+
 #if !defined(ARMA_CHECK_CONFORMANCE)
  #define ARMA_CHECK_CONFORMANCE
  //// Comment out the above line to disable conformance checks for bounds and size.
@@ -211,7 +216,7 @@


 #if defined(ARMA_DEFAULT_OSTREAM)
-  #pragma message ("WARNING: support for ARMA_DEFAULT_OSTREAM is deprecated and will be removed;")
+  #pragma message ("WARNING: option ARMA_DEFAULT_OSTREAM is deprecated and will be removed;")
  #pragma message ("WARNING: use ARMA_COUT_STREAM and ARMA_CERR_STREAM instead")
 #endif

@@ -290,7 +295,7 @@

 // for compatibility with earlier versions of Armadillo
 #if defined(ARMA_DONT_USE_CXX11_MUTEX)
-  #pragma message ("WARNING: support for ARMA_DONT_USE_CXX11_MUTEX is deprecated and will be removed;")
+  #pragma message ("WARNING: option ARMA_DONT_USE_CXX11_MUTEX is deprecated and will be removed;")
  #pragma message ("WARNING: use ARMA_DONT_USE_STD_MUTEX instead")
  #undef ARMA_USE_STD_MUTEX
 #endif
@@ -303,6 +308,18 @@
  #undef ARMA_64BIT_WORD
 #endif

+#if (defined(ARMA_BLAS_LONG_LONG) && defined(ARMA_USE_WRAPPER))
+  #pragma message ("WARNING: use of ARMA_BLAS_LONG_LONG in conjunction with ARMA_USE_WRAPPER is not supported")
+#endif
+
+#if (defined(ARMA_BLAS_64BIT_INT) && defined(ARMA_USE_WRAPPER))
+  #pragma message ("WARNING: use of ARMA_BLAS_64BIT_INT in conjunction with ARMA_USE_WRAPPER is not supported")
+#endif
+
+#if (defined(ARMA_SUPERLU_64BIT_INT) && defined(ARMA_USE_WRAPPER))
+  #pragma message ("WARNING: use of ARMA_SUPERLU_64BIT_INT in conjunction with ARMA_USE_WRAPPER is not supported")
+#endif
+
 // for compatibility with earlier versions of Armadillo
 #if defined(ARMA_BLAS_LONG) || defined(ARMA_BLAS_LONG_LONG)
  #undef  ARMA_BLAS_64BIT_INT
@@ -321,6 +338,10 @@
  #undef ARMA_OPTIMISE_INVEXPR
 #endif

+#if defined(ARMA_DONT_OPTIMISE_POWEXPR)
+  #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")
@@ -366,12 +387,12 @@
  #undef ARMA_CRIPPLED_LAPACK
 #endif

-// WARNING: option ARMA_IGNORE_DEPRECATED_MARKER will be removed
-// WARNING: option ARMA_CRIPPLED_LAPACK          will be removed
+// WARNING: option ARMA_IGNORE_DEPRECATED_MARKER is deprecated and will be removed
+// WARNING: option ARMA_CRIPPLED_LAPACK          is deprecated and will be removed

 #if defined(ARMA_CRIPPLED_LAPACK)
  #if (!defined(ARMA_IGNORE_DEPRECATED_MARKER))
-    #pragma message ("option ARMA_CRIPPLED_LAPACK is deprecated and will be removed")
+    #pragma message ("WARNING: option ARMA_CRIPPLED_LAPACK is deprecated and will be removed")
  #endif
 #endif

@@ -1441,6 +1441,7 @@ inline void arma_debug_check(bool state, const char* msg)  { arma_conform_check(
        out << "\n@ arma_config::optimise_band    = " << arma_config::optimise_band;
        out << "\n@ arma_config::optimise_sym     = " << arma_config::optimise_sym;
        out << "\n@ arma_config::optimise_invexpr = " << arma_config::optimise_invexpr;
+        out << "\n@ arma_config::optimise_powexpr = " << arma_config::optimise_powexpr;
        out << "\n@ arma_config::check_conform    = " << arma_config::check_conform;
        out << "\n@ arma_config::check_nonfinite  = " << arma_config::check_nonfinite;
        out << "\n@ arma_config::fast_math        = " << arma_config::fast_math;
@@ -293,6 +293,11 @@
  #define arma_checon checon
  #define arma_zhecon zhecon
  
+  #define arma_sgebal sgebal
+  #define arma_dgebal dgebal
+  #define arma_cgebal cgebal
+  #define arma_zgebal zgebal
+  
 #else
  
  #define arma_sgetrf SGETRF
@@ -553,6 +558,11 @@
  #define arma_checon CHECON
  #define arma_zhecon ZHECON
  
+  #define arma_sgebal SGEBAL
+  #define arma_dgebal DGEBAL
+  #define arma_cgebal CGEBAL
+  #define arma_zgebal ZGEBAL
+  
 #endif


@@ -926,6 +936,12 @@ extern "C"
  void arma_fortran(arma_checon)(const char* uplo, const blas_int* n, const blas_cxf* a, const blas_int* lda, const blas_int* ipiv, const float*  anorm, float*  rcond, blas_cxf* work, blas_int* info, blas_len uplo_len) ARMA_NOEXCEPT;
  void arma_fortran(arma_zhecon)(const char* uplo, const blas_int* n, const blas_cxd* a, const blas_int* lda, const blas_int* ipiv, const double* anorm, double* rcond, blas_cxd* work, blas_int* info, blas_len uplo_len) ARMA_NOEXCEPT;
  
+  // matrix balance
+  void arma_fortran(arma_sgebal)(const char* job, const blas_int* n,    float* a, const blas_int* lda, blas_int* ilo, blas_int* ihi,  float* scale, blas_int* info, blas_len job_len) ARMA_NOEXCEPT;
+  void arma_fortran(arma_dgebal)(const char* job, const blas_int* n,   double* a, const blas_int* lda, blas_int* ilo, blas_int* ihi, double* scale, blas_int* info, blas_len job_len) ARMA_NOEXCEPT;
+  void arma_fortran(arma_cgebal)(const char* job, const blas_int* n, blas_cxf* a, const blas_int* lda, blas_int* ilo, blas_int* ihi,  float* scale, blas_int* info, blas_len job_len) ARMA_NOEXCEPT;
+  void arma_fortran(arma_zgebal)(const char* job, const blas_int* n, blas_cxd* a, const blas_int* lda, blas_int* ilo, blas_int* ihi, double* scale, blas_int* info, blas_len job_len) ARMA_NOEXCEPT;
+  
 #else
  
  // prototypes without hidden arguments
@@ -1282,6 +1298,12 @@ extern "C"
  void arma_fortran(arma_checon)(const char* uplo, const blas_int* n, const blas_cxf* a, const blas_int* lda, const blas_int* ipiv, const float*  anorm, float*  rcond, blas_cxf* work, blas_int* info) ARMA_NOEXCEPT;
  void arma_fortran(arma_zhecon)(const char* uplo, const blas_int* n, const blas_cxd* a, const blas_int* lda, const blas_int* ipiv, const double* anorm, double* rcond, blas_cxd* work, blas_int* info) ARMA_NOEXCEPT;
  
+  // matrix balance
+  void arma_fortran(arma_sgebal)(const char* job, const blas_int* n,    float* a, const blas_int* lda, blas_int* ilo, blas_int* ihi,  float* scale, blas_int* info) ARMA_NOEXCEPT;
+  void arma_fortran(arma_dgebal)(const char* job, const blas_int* n,   double* a, const blas_int* lda, blas_int* ilo, blas_int* ihi, double* scale, blas_int* info) ARMA_NOEXCEPT;
+  void arma_fortran(arma_cgebal)(const char* job, const blas_int* n, blas_cxf* a, const blas_int* lda, blas_int* ilo, blas_int* ihi,  float* scale, blas_int* info) ARMA_NOEXCEPT;
+  void arma_fortran(arma_zgebal)(const char* job, const blas_int* n, blas_cxd* a, const blas_int* lda, blas_int* ilo, blas_int* ihi, double* scale, blas_int* info) ARMA_NOEXCEPT;
+  
 #endif
 }

@@ -19,20 +19,20 @@

 extern "C"
  {
-  extern void arma_wrapper(sgssv)(superlu::superlu_options_t*, superlu::SuperMatrix*, int*, int*, superlu::SuperMatrix*, superlu::SuperMatrix*, superlu::SuperMatrix*, superlu::SuperLUStat_t*, int*);
-  extern void arma_wrapper(dgssv)(superlu::superlu_options_t*, superlu::SuperMatrix*, int*, int*, superlu::SuperMatrix*, superlu::SuperMatrix*, superlu::SuperMatrix*, superlu::SuperLUStat_t*, int*);
-  extern void arma_wrapper(cgssv)(superlu::superlu_options_t*, superlu::SuperMatrix*, int*, int*, superlu::SuperMatrix*, superlu::SuperMatrix*, superlu::SuperMatrix*, superlu::SuperLUStat_t*, int*);
-  extern void arma_wrapper(zgssv)(superlu::superlu_options_t*, superlu::SuperMatrix*, int*, int*, superlu::SuperMatrix*, superlu::SuperMatrix*, superlu::SuperMatrix*, superlu::SuperLUStat_t*, int*);
+  extern void arma_wrapper(sgssv)(superlu::superlu_options_t*, superlu::SuperMatrix*, int*, int*, superlu::SuperMatrix*, superlu::SuperMatrix*, superlu::SuperMatrix*, superlu::SuperLUStat_t*, superlu::int_t*);
+  extern void arma_wrapper(dgssv)(superlu::superlu_options_t*, superlu::SuperMatrix*, int*, int*, superlu::SuperMatrix*, superlu::SuperMatrix*, superlu::SuperMatrix*, superlu::SuperLUStat_t*, superlu::int_t*);
+  extern void arma_wrapper(cgssv)(superlu::superlu_options_t*, superlu::SuperMatrix*, int*, int*, superlu::SuperMatrix*, superlu::SuperMatrix*, superlu::SuperMatrix*, superlu::SuperLUStat_t*, superlu::int_t*);
+  extern void arma_wrapper(zgssv)(superlu::superlu_options_t*, superlu::SuperMatrix*, int*, int*, superlu::SuperMatrix*, superlu::SuperMatrix*, superlu::SuperMatrix*, superlu::SuperLUStat_t*, superlu::int_t*);
  
-  extern void arma_wrapper(sgssvx)(superlu::superlu_options_t*, superlu::SuperMatrix*, int*, int*, int*, char*,  float*,  float*, superlu::SuperMatrix*, superlu::SuperMatrix*, void*, int, superlu::SuperMatrix*, superlu::SuperMatrix*,  float*,  float*,  float*,  float*, superlu::GlobalLU_t*, superlu::mem_usage_t*, superlu::SuperLUStat_t*, int*);
-  extern void arma_wrapper(dgssvx)(superlu::superlu_options_t*, superlu::SuperMatrix*, int*, int*, int*, char*, double*, double*, superlu::SuperMatrix*, superlu::SuperMatrix*, void*, int, superlu::SuperMatrix*, superlu::SuperMatrix*, double*, double*, double*, double*, superlu::GlobalLU_t*, superlu::mem_usage_t*, superlu::SuperLUStat_t*, int*);
-  extern void arma_wrapper(cgssvx)(superlu::superlu_options_t*, superlu::SuperMatrix*, int*, int*, int*, char*,  float*,  float*, superlu::SuperMatrix*, superlu::SuperMatrix*, void*, int, superlu::SuperMatrix*, superlu::SuperMatrix*,  float*,  float*,  float*,  float*, superlu::GlobalLU_t*, superlu::mem_usage_t*, superlu::SuperLUStat_t*, int*);
-  extern void arma_wrapper(zgssvx)(superlu::superlu_options_t*, superlu::SuperMatrix*, int*, int*, int*, char*, double*, double*, superlu::SuperMatrix*, superlu::SuperMatrix*, void*, int, superlu::SuperMatrix*, superlu::SuperMatrix*, double*, double*, double*, double*, superlu::GlobalLU_t*, superlu::mem_usage_t*, superlu::SuperLUStat_t*, int*);
+  extern void arma_wrapper(sgssvx)(superlu::superlu_options_t*, superlu::SuperMatrix*, int*, int*, int*, char*,  float*,  float*, superlu::SuperMatrix*, superlu::SuperMatrix*, void*, superlu::int_t, superlu::SuperMatrix*, superlu::SuperMatrix*,  float*,  float*,  float*,  float*, superlu::GlobalLU_t*, superlu::mem_usage_t*, superlu::SuperLUStat_t*, superlu::int_t*);
+  extern void arma_wrapper(dgssvx)(superlu::superlu_options_t*, superlu::SuperMatrix*, int*, int*, int*, char*, double*, double*, superlu::SuperMatrix*, superlu::SuperMatrix*, void*, superlu::int_t, superlu::SuperMatrix*, superlu::SuperMatrix*, double*, double*, double*, double*, superlu::GlobalLU_t*, superlu::mem_usage_t*, superlu::SuperLUStat_t*, superlu::int_t*);
+  extern void arma_wrapper(cgssvx)(superlu::superlu_options_t*, superlu::SuperMatrix*, int*, int*, int*, char*,  float*,  float*, superlu::SuperMatrix*, superlu::SuperMatrix*, void*, superlu::int_t, superlu::SuperMatrix*, superlu::SuperMatrix*,  float*,  float*,  float*,  float*, superlu::GlobalLU_t*, superlu::mem_usage_t*, superlu::SuperLUStat_t*, superlu::int_t*);
+  extern void arma_wrapper(zgssvx)(superlu::superlu_options_t*, superlu::SuperMatrix*, int*, int*, int*, char*, double*, double*, superlu::SuperMatrix*, superlu::SuperMatrix*, void*, superlu::int_t, superlu::SuperMatrix*, superlu::SuperMatrix*, double*, double*, double*, double*, superlu::GlobalLU_t*, superlu::mem_usage_t*, superlu::SuperLUStat_t*, superlu::int_t*);
  
-  extern void arma_wrapper(sgstrf)(superlu::superlu_options_t*, superlu::SuperMatrix*, int, int, int*, void*, int, int*, int*, superlu::SuperMatrix*, superlu::SuperMatrix*, superlu::GlobalLU_t*, superlu::SuperLUStat_t*, int*);
-  extern void arma_wrapper(dgstrf)(superlu::superlu_options_t*, superlu::SuperMatrix*, int, int, int*, void*, int, int*, int*, superlu::SuperMatrix*, superlu::SuperMatrix*, superlu::GlobalLU_t*, superlu::SuperLUStat_t*, int*);
-  extern void arma_wrapper(cgstrf)(superlu::superlu_options_t*, superlu::SuperMatrix*, int, int, int*, void*, int, int*, int*, superlu::SuperMatrix*, superlu::SuperMatrix*, superlu::GlobalLU_t*, superlu::SuperLUStat_t*, int*);
-  extern void arma_wrapper(zgstrf)(superlu::superlu_options_t*, superlu::SuperMatrix*, int, int, int*, void*, int, int*, int*, superlu::SuperMatrix*, superlu::SuperMatrix*, superlu::GlobalLU_t*, superlu::SuperLUStat_t*, int*);
+  extern void arma_wrapper(sgstrf)(superlu::superlu_options_t*, superlu::SuperMatrix*, int, int, int*, void*, superlu::int_t, int*, int*, superlu::SuperMatrix*, superlu::SuperMatrix*, superlu::GlobalLU_t*, superlu::SuperLUStat_t*, superlu::int_t*);
+  extern void arma_wrapper(dgstrf)(superlu::superlu_options_t*, superlu::SuperMatrix*, int, int, int*, void*, superlu::int_t, int*, int*, superlu::SuperMatrix*, superlu::SuperMatrix*, superlu::GlobalLU_t*, superlu::SuperLUStat_t*, superlu::int_t*);
+  extern void arma_wrapper(cgstrf)(superlu::superlu_options_t*, superlu::SuperMatrix*, int, int, int*, void*, superlu::int_t, int*, int*, superlu::SuperMatrix*, superlu::SuperMatrix*, superlu::GlobalLU_t*, superlu::SuperLUStat_t*, superlu::int_t*);
+  extern void arma_wrapper(zgstrf)(superlu::superlu_options_t*, superlu::SuperMatrix*, int, int, int*, void*, superlu::int_t, int*, int*, superlu::SuperMatrix*, superlu::SuperMatrix*, superlu::GlobalLU_t*, superlu::SuperLUStat_t*, superlu::int_t*);
  
  extern void arma_wrapper(sgstrs)(superlu::trans_t, superlu::SuperMatrix*, superlu::SuperMatrix*, int*, int*, superlu::SuperMatrix*, superlu::SuperLUStat_t*, int*);
  extern void arma_wrapper(dgstrs)(superlu::trans_t, superlu::SuperMatrix*, superlu::SuperMatrix*, int*, int*, superlu::SuperMatrix*, superlu::SuperLUStat_t*, int*);
@@ -995,9 +995,13 @@ diagview<eT>::randu()
  
  const uword local_n_elem = n_elem;
  
+  Col<eT> tmp(local_n_elem, arma_nozeros_indicator());
+  
+  tmp.randu();
+  
  for(uword ii=0; ii < local_n_elem; ++ii)
    {
-    x.at(ii+row_offset, ii+col_offset) = eT(arma_rng::randu<eT>());
+    x.at(ii+row_offset, ii+col_offset) = tmp[ii];
    }
  }

@@ -1014,9 +1018,13 @@ diagview<eT>::randn()
  
  const uword local_n_elem = n_elem;
  
+  Col<eT> tmp(local_n_elem, arma_nozeros_indicator());
+  
+  tmp.randn();
+  
  for(uword ii=0; ii < local_n_elem; ++ii)
    {
-    x.at(ii+row_offset, ii+col_offset) = eT(arma_rng::randn<eT>());
+    x.at(ii+row_offset, ii+col_offset) = tmp[ii];
    }
  }

@@ -2120,7 +2120,12 @@ field<oT>::init(const field<oT>& x)
  {
  arma_debug_sigprint();
  
-  if(this == &x)  { return; }
+  if(this == &x)
+    {
+    arma_debug_print("field::init(): copy omitted");
+    
+    return;
+    }
  
  field& t = (*this);
  
@@ -272,6 +272,99 @@ accu(const T1& X)



+template<typename T1, typename functor>
+inline
+typename T1::elem_type
+accu_op_omit_helper(const Proxy<T1>& P, functor is_omitted)
+  {
+  arma_debug_sigprint();
+  
+  typedef typename T1::elem_type eT;
+  
+  constexpr eT eT_zero = eT(0);
+  
+  eT acc = eT(0);
+  
+  if(Proxy<T1>::use_at)
+    {
+    const uword n_rows = P.get_n_rows();
+    const uword n_cols = P.get_n_cols();
+    
+    for(uword c=0; c < n_cols; ++c)
+    for(uword r=0; r < n_rows; ++r)
+      {
+      const eT val = P.at(r,c);
+      
+      acc += is_omitted(val) ? eT_zero : val;
+      }
+    }
+  else
+    {
+    typename Proxy<T1>::ea_type Pea = P.get_ea();
+    
+    const uword n_elem = P.get_n_elem();
+    
+    eT val1 = eT(0);
+    eT val2 = eT(0);
+    
+    uword i,j;
+    for(i=0, j=1; j < n_elem; i+=2, j+=2)
+      {
+      const eT tmp_i = Pea[i];
+      const eT tmp_j = Pea[j];
+      
+      val1 += is_omitted(tmp_i) ? eT_zero : tmp_i;
+      val2 += is_omitted(tmp_j) ? eT_zero : tmp_j;
+      }
+    
+    if(i < n_elem)
+      {
+      const eT tmp_i = Pea[i];
+      
+      val1 += is_omitted(tmp_i) ? eT_zero : tmp_i;
+      }
+    
+    acc = val1 + val2;
+    }
+  
+  return acc;
+  }
+
+
+
+template<typename T1>
+arma_warn_unused
+inline
+typename T1::elem_type
+accu(const Op<T1, op_omit>& in)
+  {
+  arma_debug_sigprint();
+  
+  typedef typename T1::elem_type eT;
+  
+  const uword omit_mode = in.aux_uword_a;
+  
+  if(arma_config::fast_math_warn)
+    {
+    if(omit_mode == 1)  { arma_warn(1, "omit_nan(): detection of NaN is not reliable in fast math mode"); }
+    if(omit_mode == 2)  { arma_warn(1, "omit_nonfinite(): detection of non-finite values is not reliable in fast math mode"); }
+    }
+  
+  auto is_omitted_1 = [](const eT& x) -> bool  { return arma_isnan(x);       };
+  auto is_omitted_2 = [](const eT& x) -> bool  { return arma_isnonfinite(x); };
+  
+  const Proxy<T1> P(in.m);
+  
+  eT acc = eT(0);
+  
+  if(omit_mode == 1)  { acc = accu_op_omit_helper(P, is_omitted_1); }
+  if(omit_mode == 2)  { acc = accu_op_omit_helper(P, is_omitted_2); }
+  
+  return acc;
+  }
+
+
+
 template<typename T1>
 arma_warn_unused
 inline
@@ -314,14 +407,14 @@ accu(const eOp<T1,eop_pow>& expr)
  
  typedef eOp<T1,eop_pow> expr_type;
  
-  if(expr.aux == eT(2))
+  if(arma_config::optimise_powexpr && (expr.aux == eT(2)))
    {
    typedef eOp<T1,eop_square> modified_expr_type;
    
    return accu( reinterpret_cast< const modified_expr_type& >(expr) );
    }
  
-  if((expr.aux == eT(0.5)) && is_non_integral<eT>::value)
+  if(arma_config::optimise_powexpr && (expr.aux == eT(0.5)) && is_non_integral<eT>::value)
    {
    typedef eOp<T1,eop_sqrt> modified_expr_type;
    
@@ -955,14 +1048,14 @@ accu(const eOpCube<T1,eop_pow>& expr)
  
  typedef eOpCube<T1,eop_pow> expr_type;
  
-  if(expr.aux == eT(2))
+  if(arma_config::optimise_powexpr && (expr.aux == eT(2)))
    {
    typedef eOpCube<T1,eop_square> modified_expr_type;
    
    return accu( reinterpret_cast< const modified_expr_type& >(expr) );
    }
  
-  if((expr.aux == eT(0.5)) && is_non_integral<eT>::value)
+  if(arma_config::optimise_powexpr && (expr.aux == eT(0.5)) && is_non_integral<eT>::value)
    {
    typedef eOpCube<T1,eop_sqrt> modified_expr_type;
    
@@ -1005,6 +1098,101 @@ accu(const eGlueCube<T1,T2,eglue_schur>& expr)



+template<typename T1, typename functor>
+inline
+typename T1::elem_type
+accu_cube_omit_helper(const ProxyCube<T1>& P, functor is_omitted)
+  {
+  arma_debug_sigprint();
+  
+  typedef typename T1::elem_type eT;
+  
+  constexpr eT eT_zero = eT(0);
+  
+  eT acc = eT(0);
+  
+  if(ProxyCube<T1>::use_at)
+    {
+    const uword n_r = P.get_n_rows();
+    const uword n_c = P.get_n_cols();
+    const uword n_s = P.get_n_slices();
+    
+    for(uword s=0; s < n_s; ++s)
+    for(uword c=0; c < n_c; ++c)
+    for(uword r=0; r < n_r; ++r)
+      {
+      const eT val = P.at(r,c,s);
+      
+      acc += is_omitted(val) ? eT_zero : val;
+      }
+    }
+  else
+    {
+    typename ProxyCube<T1>::ea_type Pea = P.get_ea();
+    
+    const uword n_elem = P.get_n_elem();
+    
+    eT val1 = eT(0);
+    eT val2 = eT(0);
+    
+    uword i,j;
+    for(i=0, j=1; j < n_elem; i+=2, j+=2)
+      {
+      const eT tmp_i = Pea[i];
+      const eT tmp_j = Pea[j];
+      
+      val1 += is_omitted(tmp_i) ? eT_zero : tmp_i;
+      val2 += is_omitted(tmp_j) ? eT_zero : tmp_j;
+      }
+    
+    if(i < n_elem)
+      {
+      const eT tmp_i = Pea[i];
+      
+      val1 += is_omitted(tmp_i) ? eT_zero : tmp_i;
+      }
+    
+    acc = val1 + val2;
+    }
+  
+  return acc;
+  }
+
+
+
+template<typename T1>
+arma_warn_unused
+inline
+typename T1::elem_type
+accu(const CubeToMatOp<T1, op_omit_cube>& in)
+  {
+  arma_debug_sigprint();
+  
+  typedef typename T1::elem_type eT;
+  
+  const ProxyCube<T1> P(in.m);
+  
+  const uword omit_mode = in.aux_uword;
+  
+  if(arma_config::fast_math_warn)
+    {
+    if(omit_mode == 1)  { arma_warn(1, "omit_nan(): detection of NaN is not reliable in fast math mode"); }
+    if(omit_mode == 2)  { arma_warn(1, "omit_nonfinite(): detection of non-finite values is not reliable in fast math mode"); }
+    }
+  
+  auto is_omitted_1 = [](const eT& x) -> bool  { return arma_isnan(x);       };
+  auto is_omitted_2 = [](const eT& x) -> bool  { return arma_isnonfinite(x); };
+  
+  eT acc = eT(0);
+  
+  if(omit_mode == 1)  { acc = accu_cube_omit_helper(P, is_omitted_1); }
+  if(omit_mode == 2)  { acc = accu_cube_omit_helper(P, is_omitted_2); }
+  
+  return acc;
+  }
+
+
+
 //


@@ -1121,6 +1309,8 @@ accu(const SpGlue<T1,T2,spglue_schur>& expr)
  
  arma_conform_assert_same_size(px.get_n_rows(), px.get_n_cols(), py.get_n_rows(), py.get_n_cols(), "element-wise multiplication");
  
+  if( (px.get_n_nonzero() == 0) && (py.get_n_nonzero() == 0) )  { return eT(0); }
+  
  typedef typename SpProxy<T1>::stored_type px_Q_type;
  typedef typename SpProxy<T2>::stored_type py_Q_type;
  
@@ -1162,10 +1352,14 @@ accu(const SpGlue<T1,T2,spglue_schur>& expr)
      
      if((x_it_col < y_it_col) || ((x_it_col == y_it_col) && (x_it_row < y_it_row))) // if y is closer to the end
        {
+        acc += (*x_it) * eT(0);  // in case (*x_it) is inf or nan
+        
        ++x_it;
        }
      else // x is closer to the end
        {
+        acc += eT(0) * (*y_it);  // in case (*y_it) is inf or nan
+        
        ++y_it;
        }
      }
@@ -1193,30 +1387,6 @@ accu(const SpOp<T1, spop_type>& expr)
  
  if(is_vectorise)  { return accu(expr.m); }
  
-  if(is_same_type<spop_type, spop_square>::yes)
-    {
-    const SpProxy<T1> P(expr.m);
-    
-    const uword N = P.get_n_nonzero();
-    
-    if(N == 0)  { return eT(0); }
-    
-    if(SpProxy<T1>::use_iterator == false)
-      {
-      return op_dot::direct_dot(N, P.get_values(), P.get_values());
-      }
-    else
-      {
-      typename SpProxy<T1>::const_iterator_type it = P.begin();
-      
-      eT val = eT(0);
-      
-      for(uword i=0; i < N; ++i)  { const eT tmp = (*it); val += (tmp*tmp); ++it; }
-      
-      return val;
-      }
-    }
-  
  const SpMat<eT> tmp = expr;
  
  return accu(tmp);
@@ -1224,6 +1394,169 @@ accu(const SpOp<T1, spop_type>& expr)



+template<typename T1>
+arma_warn_unused
+inline
+typename T1::elem_type
+accu(const SpOp<T1, spop_square>& expr)
+  {
+  arma_debug_sigprint();
+  
+  typedef typename T1::elem_type eT;
+  
+  if(is_SpSubview_col<T1>::value)
+    {
+    const SpSubview_col<eT>& svcol = reinterpret_cast<const SpSubview_col<eT>&>(expr.m);
+    
+    if(svcol.n_nonzero == 0)  { return eT(0); }
+    
+    if(svcol.n_rows == svcol.m.n_rows)
+      {
+      arma_debug_print("accu(): SpSubview_col spop_square optimisation");
+      
+      const SpMat<eT>& m   = svcol.m;
+      const uword      col = svcol.aux_col1;
+      
+      const eT* ptr = &(m.values[ m.col_ptrs[col] ]);
+      
+      return op_dot::direct_dot(svcol.n_nonzero, ptr, ptr);
+      }
+    }
+  
+  const SpProxy<T1> P(expr.m);
+  
+  const uword N = P.get_n_nonzero();
+  
+  if(N == 0)  { return eT(0); }
+  
+  if(SpProxy<T1>::use_iterator == false)
+    {
+    return op_dot::direct_dot(N, P.get_values(), P.get_values());
+    }
+  else
+    {
+    typename SpProxy<T1>::const_iterator_type it = P.begin();
+    
+    eT acc = eT(0);
+    
+    for(uword i=0; i < N; ++i)  { const eT tmp = (*it); acc += (tmp*tmp); ++it; }
+    
+    return acc;
+    }
+  }
+
+
+
+template<typename T1, typename functor>
+inline
+typename T1::elem_type
+accu_spop_omit_helper(const T1& expr, functor is_omitted)
+  {
+  arma_debug_sigprint();
+  
+  typedef typename T1::elem_type eT;
+  
+  constexpr eT eT_zero = eT(0);
+  
+  if(is_SpSubview_col<T1>::value)
+    {
+    const SpSubview_col<eT>& svcol = reinterpret_cast<const SpSubview_col<eT>&>(expr);
+    
+    if(svcol.n_nonzero == 0)  { return eT(0); }
+    
+    if(svcol.n_rows == svcol.m.n_rows)
+      {
+      arma_debug_print("accu_spop_omit_helper(): SpSubview_col optimisation");
+      
+      const SpMat<eT>& m   = svcol.m;
+      const uword      col = svcol.aux_col1;
+      
+      const eT* vals = &(m.values[ m.col_ptrs[col] ]);
+      
+      const uword N = svcol.n_nonzero;
+      
+      eT acc = eT(0);
+      
+      for(uword i=0; i < N; ++i)
+        {
+        const eT tmp = vals[i];
+        
+        acc += is_omitted(tmp) ? eT_zero : tmp;
+        }
+      
+      return acc;
+      }
+    }
+  
+  const SpProxy<T1> P(expr);
+  
+  const uword N = P.get_n_nonzero();
+  
+  if(N == 0)  { return eT(0); }
+  
+  eT acc = eT(0);
+  
+  if(SpProxy<T1>::use_iterator == false)
+    {
+    const eT* vals = P.get_values();
+    
+    for(uword i=0; i < N; ++i)
+      {
+      const eT tmp = vals[i];
+      
+      acc += is_omitted(tmp) ? eT_zero : tmp;
+      }
+    }
+  else
+    {
+    typename SpProxy<T1>::const_iterator_type it = P.begin();
+    
+    for(uword i=0; i < N; ++i)
+      {
+      const eT tmp = (*it);
+      
+      acc += is_omitted(tmp) ? eT_zero : tmp;
+      
+      ++it;
+      }
+    }
+  
+  return acc;
+  }
+
+
+
+template<typename T1>
+arma_warn_unused
+inline
+typename T1::elem_type
+accu(const SpOp<T1, spop_omit>& expr)
+  {
+  arma_debug_sigprint();
+  
+  typedef typename T1::elem_type eT;
+  
+  const uword omit_mode = expr.aux_uword_a;
+  
+  if(arma_config::fast_math_warn)
+    {
+    if(omit_mode == 1)  { arma_warn(1, "omit_nan(): detection of NaN is not reliable in fast math mode"); }
+    if(omit_mode == 2)  { arma_warn(1, "omit_nonfinite(): detection of non-finite values is not reliable in fast math mode"); }
+    }
+  
+  auto is_omitted_1 = [](const eT& x) -> bool { return arma_isnan(x);       };
+  auto is_omitted_2 = [](const eT& x) -> bool { return arma_isnonfinite(x); };
+  
+  eT acc = eT(0);
+  
+  if(omit_mode == 1)  { acc = accu_spop_omit_helper(expr.m, is_omitted_1); }
+  if(omit_mode == 2)  { acc = accu_spop_omit_helper(expr.m, is_omitted_2); }
+  
+  return acc;
+  }
+
+
+
 template<typename T1, typename spop_type>
 arma_warn_unused
 inline
@@ -0,0 +1,135 @@
+// SPDX-License-Identifier: Apache-2.0
+// 
+// Copyright 2008-2016 Conrad Sanderson (http://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
+// http://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_balance
+//! @{
+
+
+
+template<typename T1>
+inline
+typename enable_if2< is_supported_blas_type<typename T1::elem_type>::value, bool >::result
+balance(Col<typename T1::pod_type>& S, Col<uword>& P, Mat<typename T1::elem_type>& B, const Base<typename T1::elem_type,T1>& A, const char* method = "both")
+  {
+  arma_debug_sigprint();
+  
+  arma_conform_check( (void_ptr(&S)  == void_ptr(&B)), "eig_gen(): parameter 'S' is an alias of parameter 'B'" );
+  
+  const char sig = (method != nullptr) ? method[0] : char(0);
+  
+  if( (sig != 'b') && (sig != 's') && (sig != 'p') )  { arma_stop_logic_error("balance(): unsupported method"); }
+  
+  const bool do_scale = (sig == 'b') || (sig == 's');
+  const bool do_perm  = (sig == 'b') || (sig == 'p');
+  
+  const bool calc_SP = true;
+  
+  B = A.get_ref();
+  
+  if(B.is_square() == false)
+    {
+    B.soft_reset();
+    
+    arma_stop_logic_error("balance(): given matrix must be square sized");
+    
+    return false;
+    }
+  
+  const bool status = auxlib::balance(S, P, B, calc_SP, do_scale, do_perm);
+  
+  if(status == false)
+    {
+    S.soft_reset();
+    P.soft_reset();
+    B.soft_reset();
+    
+    arma_warn(3, "balance(): transformation failed");
+    }
+  
+  return status;
+  }
+
+
+
+template<typename T1>
+inline
+typename enable_if2< is_supported_blas_type<typename T1::elem_type>::value, bool >::result
+balance(Mat<typename T1::elem_type>& B, const Base<typename T1::elem_type,T1>& A, const char* method = "both")
+  {
+  arma_debug_sigprint();
+  
+  typedef typename T1::pod_type T;
+  
+  const char sig = (method != nullptr) ? method[0] : char(0);
+  
+  if( (sig != 'b') && (sig != 's') && (sig != 'p') )  { arma_stop_logic_error("balance(): unsupported method"); }
+  
+  const bool do_scale = (sig == 'b') || (sig == 's');
+  const bool do_perm  = (sig == 'b') || (sig == 'p');
+  
+  const bool calc_SP = false;
+  
+  B = A.get_ref();
+  
+  if(B.is_square() == false)
+    {
+    B.soft_reset();
+    
+    arma_stop_logic_error("balance(): given matrix must be square sized");
+    
+    return false;
+    }
+  
+  Col<T>     S;
+  Col<uword> P;
+  
+  const bool status = auxlib::balance(S, P, B, calc_SP, do_scale, do_perm);
+  
+  if(status == false)
+    {
+    B.soft_reset();
+    
+    arma_warn(3, "balance(): transformation failed");
+    }
+  
+  return status;
+  }
+
+
+
+template<typename T1>
+inline
+typename enable_if2< is_supported_blas_type<typename T1::elem_type>::value, Mat<typename T1::elem_type> >::result
+balance(const Base<typename T1::elem_type,T1>& A, const char* method = "both")
+  {
+  arma_debug_sigprint();
+  
+  typedef typename T1::elem_type eT;
+  
+  Mat<eT> B;
+  
+  const bool status = balance(B, A, method);
+  
+  if(status == false)  { arma_stop_runtime_error("balance(): transformation failed"); }
+  
+  return B;
+  }
+
+
+
+//! @}
@@ -29,19 +29,19 @@ class conv_to
  public:
  
  template<typename in_eT>
-  arma_frown("use as_scalar() instead") inline static out_eT from(const in_eT& in, const typename arma_scalar_only<in_eT>::result* junk = nullptr);
+  arma_frown("replace conv_to<...>::from(X) with as_scalar(X)") inline static out_eT from(const in_eT& in, const typename arma_scalar_only<in_eT>::result* junk = nullptr);
  
  template<typename in_eT, typename T1>
-  arma_frown("use as_scalar() instead") inline static out_eT from(const Base<in_eT, T1>& in, const typename arma_not_cx<in_eT>::result* junk = nullptr);
+  arma_frown("replace conv_to<...>::from(X) with as_scalar(X)") inline static out_eT from(const Base<in_eT, T1>& in, const typename arma_not_cx<in_eT>::result* junk = nullptr);
  
  template<typename in_eT, typename T1>
-  arma_frown("use as_scalar() instead") inline static out_eT from(const Base<in_eT, T1>& in, const typename arma_cx_only<in_eT>::result* junk = nullptr);
+  arma_frown("replace conv_to<...>::from(X) with as_scalar(X)") inline static out_eT from(const Base<in_eT, T1>& in, const typename arma_cx_only<in_eT>::result* junk = nullptr);
  
  template<typename in_eT, typename T1>
-  arma_frown("use as_scalar() instead") inline static out_eT from(const BaseCube<in_eT, T1>& in, const typename arma_not_cx<in_eT>::result* junk = nullptr);
+  arma_frown("replace conv_to<...>::from(X) with as_scalar(X)") inline static out_eT from(const BaseCube<in_eT, T1>& in, const typename arma_not_cx<in_eT>::result* junk = nullptr);
  
  template<typename in_eT, typename T1>
-  arma_frown("use as_scalar() instead") inline static out_eT from(const BaseCube<in_eT, T1>& in, const typename arma_cx_only<in_eT>::result* junk = nullptr);
+  arma_frown("replace conv_to<...>::from(X) with as_scalar(X)") inline static out_eT from(const BaseCube<in_eT, T1>& in, const typename arma_cx_only<in_eT>::result* junk = nullptr);
  };


@@ -545,6 +545,45 @@ arg(const T1& X)



+template<typename T1>
+arma_warn_unused
+inline
+typename enable_if2< is_arma_type<T1>::value, const mtOp<typename T1::elem_type, T1, op_replace> >::result
+replace(const T1& X, typename T1::elem_type old_val, typename T1::elem_type new_val)
+  {
+  arma_debug_sigprint();
+  
+  return mtOp<typename T1::elem_type, T1, op_replace>(mtOp_dual_aux_indicator(), X, old_val, new_val);
+  }
+
+
+
+template<typename T1>
+arma_warn_unused
+inline
+const mtOpCube<typename T1::elem_type, T1, op_replace>
+replace(const BaseCube<typename T1::elem_type,T1>& X, typename T1::elem_type old_val, typename T1::elem_type new_val)
+  {
+  arma_debug_sigprint();
+  
+  return mtOpCube<typename T1::elem_type, T1, op_replace>(mtOpCube_dual_aux_indicator(), X.get_ref(), old_val, new_val);
+  }
+
+
+
+template<typename T1>
+arma_warn_unused
+inline
+typename enable_if2< is_arma_sparse_type<T1>::value, const mtSpOp<typename T1::elem_type, T1, spop_replace> >::result
+replace(const T1& X, typename T1::elem_type old_val, typename T1::elem_type new_val)
+  {
+  arma_debug_sigprint();
+  
+  return mtSpOp<typename T1::elem_type, T1, spop_replace>(mtSpOp_dual_aux_indicator(), X, old_val, new_val);
+  }
+
+
+
 //
 // square

@@ -269,6 +269,23 @@ find_nan(const T1& X)



+template<typename T1>
+arma_warn_unused
+inline
+typename enable_if2
+  <
+  is_arma_type<T1>::value,
+  const mtOp<uword, T1, op_find_nonnan>
+  >::result
+find_nonnan(const T1& X)
+  {
+  arma_debug_sigprint();
+
+  return mtOp<uword, T1, op_find_nonnan>(X);
+  }
+
+
+
 //


@@ -330,6 +347,25 @@ find_nan(const BaseCube<typename T1::elem_type,T1>& X)



+template<typename T1>
+arma_warn_unused
+inline
+uvec
+find_nonnan(const BaseCube<typename T1::elem_type,T1>& X)
+  {
+  arma_debug_sigprint();
+
+  typedef typename T1::elem_type eT;
+
+  const unwrap_cube<T1> tmp(X.get_ref());
+
+  const Mat<eT> R( const_cast< eT* >(tmp.M.memptr()), tmp.M.n_elem, 1, false );
+
+  return find_nonnan(R);
+  }
+
+
+
 //


@@ -401,7 +437,7 @@ find_nonfinite(const SpBase<typename T1::elem_type,T1>& X)
  
  for(uword i=0; i<n_nz; ++i)
    {
-    if(arma_isfinite(*it) == false)
+    if(arma_isnonfinite(*it))
      {
      const uword index = it.row() + it.col()*n_rows;
      
@@ -466,4 +502,48 @@ find_nan(const SpBase<typename T1::elem_type,T1>& X)



+template<typename T1>
+arma_warn_unused
+inline
+Col<uword>
+find_nonnan(const SpBase<typename T1::elem_type,T1>& X)
+  {
+  arma_debug_sigprint();
+
+  const SpProxy<T1> P(X.get_ref());
+
+  const uword n_rows = P.get_n_rows();
+  const uword n_nz   = P.get_n_nonzero();
+
+  Mat<uword> tmp(n_nz, 1, arma_nozeros_indicator());
+
+  uword* tmp_mem = tmp.memptr();
+
+  typename SpProxy<T1>::const_iterator_type it = P.begin();
+
+  uword count = 0;
+
+  for(uword i=0; i<n_nz; ++i)
+    {
+    if(arma_isnan(*it) == false)
+      {
+      const uword index = it.row() + it.col()*n_rows;
+
+      tmp_mem[count] = index;
+
+      ++count;
+      }
+
+    ++it;
+    }
+
+  Col<uword> out;
+
+  if(count > 0)  { out.steal_mem_col(tmp, count); }
+
+  return out;
+  }
+
+
+
 //! @}
@@ -149,8 +149,8 @@ logspace(const double A, const double B, const uword N = 50u)

 //! kept for compatibility with old user code
 template<typename eT>
-arma_warn_unused
-arma_inline
+arma_frown("change arma::is_finite(val) to std::isfinite(val)")
+inline
 bool
 is_finite(const eT x, const typename arma_scalar_only<eT>::result* junk = nullptr)
  {
@@ -163,7 +163,7 @@ is_finite(const eT x, const typename arma_scalar_only<eT>::result* junk = nullpt

 //! kept for compatibility with old user code
 template<typename T1>
-arma_warn_unused
+arma_frown("change arma::is_finite(X) to X.is_finite()")
 inline
 bool
 is_finite(const Base<typename T1::elem_type,T1>& X)
@@ -177,7 +177,7 @@ is_finite(const Base<typename T1::elem_type,T1>& X)

 //! kept for compatibility with old user code
 template<typename T1>
-arma_warn_unused
+arma_frown("change arma::is_finite(X) to X.is_finite()")
 inline
 bool
 is_finite(const SpBase<typename T1::elem_type,T1>& X)
@@ -191,7 +191,7 @@ is_finite(const SpBase<typename T1::elem_type,T1>& X)

 //! kept for compatibility with old user code
 template<typename T1>
-arma_warn_unused
+arma_frown("change arma::is_finite(X) to X.is_finite()")
 inline
 bool
 is_finite(const BaseCube<typename T1::elem_type,T1>& X)
@@ -557,7 +557,7 @@ namespace priv
    
    const eT negdelta = log_b - log_a;
    
-    if( (negdelta < Datum<eT>::log_min) || (arma_isfinite(negdelta) == false) )
+    if( (negdelta < Datum<eT>::log_min) || arma_isnonfinite(negdelta) )
      {
      return log_a;
      }
@@ -0,0 +1,104 @@
+// SPDX-License-Identifier: Apache-2.0
+// 
+// Copyright 2008-2016 Conrad Sanderson (http://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
+// http://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_omit
+//! @{
+
+
+
+template<typename T1>
+arma_warn_unused
+inline
+typename enable_if2< is_arma_type<T1>::value, const Op<T1, op_omit> >::result
+omit_nan(const T1& X)
+  {
+  arma_debug_sigprint();
+  
+  return Op<T1, op_omit>(X, 1, 0);
+  }
+
+
+
+template<typename T1>
+arma_warn_unused
+inline
+typename enable_if2< is_arma_type<T1>::value, const Op<T1, op_omit> >::result
+omit_nonfinite(const T1& X)
+  {
+  arma_debug_sigprint();
+  
+  return Op<T1, op_omit>(X, 2, 0);
+  }
+
+
+
+template<typename T1>
+arma_warn_unused
+inline
+CubeToMatOp<T1, op_omit_cube>
+omit_nan(const BaseCube<typename T1::elem_type, T1>& X)
+  {
+  arma_debug_sigprint();
+  
+  return CubeToMatOp<T1, op_omit_cube>(X.get_ref(), 1);
+  }
+
+
+
+template<typename T1>
+arma_warn_unused
+inline
+CubeToMatOp<T1, op_omit_cube>
+omit_nonfinite(const BaseCube<typename T1::elem_type, T1>& X)
+  {
+  arma_debug_sigprint();
+  
+  return CubeToMatOp<T1, op_omit_cube>(X.get_ref(), 2);
+  }
+
+
+
+template<typename T1>
+arma_warn_unused
+inline
+typename
+enable_if2< is_arma_sparse_type<T1>::value, const SpOp<T1, spop_omit> >::result
+omit_nan(const T1& X)
+  {
+  arma_debug_sigprint();
+  
+  return SpOp<T1, spop_omit>(X, 1, 0);
+  }
+
+
+
+template<typename T1>
+arma_warn_unused
+inline
+typename
+enable_if2< is_arma_sparse_type<T1>::value, const SpOp<T1, spop_omit> >::result
+omit_nonfinite(const T1& X)
+  {
+  arma_debug_sigprint();
+  
+  return SpOp<T1, spop_omit>(X, 2, 0);
+  }
+
+
+
+//! @}
@@ -44,7 +44,7 @@ pow


 template<typename parent, unsigned int mode, typename T2>
-arma_deprecated
+arma_frown("refactor your code to use pow() in conjunction with repmat()")
 inline
 Mat<typename parent::elem_type>
 pow
@@ -121,7 +121,7 @@ pow


 template<typename parent, unsigned int mode, typename T2>
-arma_deprecated
+arma_frown("refactor your code to use pow() in conjunction with repmat()")
 inline
 typename
 enable_if2
@@ -37,9 +37,18 @@ internal_randperm_helper(obj_type& x, const uword N, const uword N_keep)
  
  std::vector<packet> packet_vec(N);
  
+  podarray<int> tmp(N);
+  
+  int* tmp_mem = tmp.memptr();
+  
+  const int a = 0;
+  const int b = arma_rng::randi<int>::max_val();
+  
+  arma_rng::randi<int>::fill(tmp_mem, N, a, b);
+  
  for(uword i=0; i < N; ++i)
    {
-    packet_vec[i].val   = int(arma_rng::randi<int>());
+    packet_vec[i].val   = tmp_mem[i];
    packet_vec[i].index = i;
    }
  
@@ -21,6 +21,9 @@
 //! @{


+// TODO: deprecate shift() in favour of circshift()
+
+
 template<typename T1>
 arma_warn_unused
 arma_inline
@@ -67,10 +67,11 @@ sort_index



+// DO NOT USE: kept only for compatibility with old user code
 template<typename T1>
 arma_warn_unused
 arma_inline
-const mtOp<uword,T1,op_stable_sort_index>
+const mtOp<uword,T1,op_sort_index>
 stable_sort_index
  (
  const Base<typename T1::elem_type,T1>& X
@@ -78,11 +79,12 @@ stable_sort_index
  {
  arma_debug_sigprint();
  
-  return mtOp<uword,T1,op_stable_sort_index>(X.get_ref(), uword(0), uword(0));
+  return mtOp<uword,T1,op_sort_index>(X.get_ref(), uword(0), uword(0));
  }



+// DO NOT USE: kept only for compatibility with old user code
 template<typename T1, typename T2>
 arma_warn_unused
 inline
@@ -90,7 +92,7 @@ typename
 enable_if2
  <
  ( (is_arma_type<T1>::value) && (is_same_type<T2, char>::value) ),
-  const mtOp<uword,T1,op_stable_sort_index>
+  const mtOp<uword,T1,op_sort_index>
  >::result
 stable_sort_index
  (
@@ -104,7 +106,7 @@ stable_sort_index
  
  arma_conform_check( ((sig != 'a') && (sig != 'd')), "stable_sort_index(): unknown sort direction" );
  
-  return mtOp<uword,T1,op_stable_sort_index>(X, ((sig == 'a') ? uword(0) : uword(1)), uword(0));
+  return mtOp<uword,T1,op_sort_index>(X, ((sig == 'a') ? uword(0) : uword(1)), uword(0));
  }


@@ -24,7 +24,7 @@
 template<typename T1, typename T2, typename T3>
 inline
 bool
-syl
+sylvester
  (
        Mat <typename T1::elem_type>   & out,
  const Base<typename T1::elem_type,T1>& in_A,
@@ -46,12 +46,12 @@ syl
  const Mat<eT>& B = tmp_B.M;
  const Mat<eT>& C = tmp_C.M;
  
-  const bool status = auxlib::syl(out, A, B, C);
+  const bool status = auxlib::sylvester(out, A, B, C);
  
  if(status == false)
    {
    out.soft_reset();
-    arma_warn(3, "syl(): solution not found");
+    arma_warn(3, "sylvester(): solution not found");
    }
  
  return status;
@@ -59,10 +59,12 @@ syl



+// kept for compatibility with old user code
 template<typename T1, typename T2, typename T3>
+arma_frown("use sylvester() instead")
 inline
 bool
-sylvester
+syl
  (
        Mat <typename T1::elem_type>   & out,
  const Base<typename T1::elem_type,T1>& in_A,
@@ -72,7 +74,8 @@ sylvester
  )
  {
  arma_ignore(junk);
-  return syl(out, in_A, in_B, in_C);
+  
+  return sylvester(out, in_A, in_B, in_C);
  }


@@ -81,7 +84,7 @@ template<typename T1, typename T2, typename T3>
 arma_warn_unused
 inline
 Mat<typename T1::elem_type>
-syl
+sylvester
  (
  const Base<typename T1::elem_type,T1>& in_A,
  const Base<typename T1::elem_type,T2>& in_B,
@@ -104,12 +107,12 @@ syl
  
  Mat<eT> out;
  
-  const bool status = auxlib::syl(out, A, B, C);
+  const bool status = auxlib::sylvester(out, A, B, C);
  
  if(status == false)
    {
    out.soft_reset();
-    arma_stop_runtime_error("syl(): solution not found");
+    arma_stop_runtime_error("sylvester(): solution not found");
    }
  
  return out;
@@ -117,11 +120,12 @@ syl



+// kept for compatibility with old user code
 template<typename T1, typename T2, typename T3>
-arma_warn_unused
+arma_frown("use sylvester() instead")
 inline
 Mat<typename T1::elem_type>
-sylvester
+syl
  (
  const Base<typename T1::elem_type,T1>& in_A,
  const Base<typename T1::elem_type,T2>& in_B,
@@ -130,7 +134,8 @@ sylvester
  )
  {
  arma_ignore(junk);
-  return syl(in_A, in_B, in_C);
+  
+  return sylvester(in_A, in_B, in_C);
  }


@@ -139,13 +139,13 @@ glue_mvnrnd::apply_noalias(Mat<eT>& out, const Mat<eT>& M, const Mat<eT>& C, con
    
    const eT tol = eT(-100) * Datum<eT>::eps * norm(C, "fro");
    
-    if(arma_isfinite(tol) == false)  { return false; }
+    if(arma_isnonfinite(tol))  { return false; }
    
    for(uword i=0; i<eigval_n_elem; ++i)
      {
      const eT val = eigval_mem[i];
      
-      if( (val < tol) || (arma_isfinite(val) == false) )  { return false; }
+      if( (val < tol) || arma_isnonfinite(val) )  { return false; }
      }
    
    for(uword i=0; i<eigval_n_elem; ++i)  { if(eigval_mem[i] < eT(0))  { eigval_mem[i] = eT(0); } }
@@ -104,6 +104,40 @@ glue_times_redirect2_helper<true>::apply(Mat<typename T1::elem_type>& out, const
    
    const strip_inv<T1> A_strip(X.A);
    
+    typedef typename strip_inv<T1>::stored_type T1_stripped;
+    
+    if( (is_cx<eT>::no) && (strip_inv<T1>::do_inv_gen) && (is_Mat<T1_stripped>::value) && (is_Mat<T2>::value) )
+      {
+      const unwrap<T1_stripped> UA(A_strip.M);
+      const unwrap<T2         > UB(X.B);
+      
+      const typename unwrap<T1_stripped>::stored_type& A = UA.M;
+      const typename unwrap<T2         >::stored_type& B = UB.M;
+      
+      const uword N = A.n_rows;
+      
+      if( (N > 0) && (N <= uword(3)) && (N == A.n_cols) && (N == B.n_rows) && (void_ptr(&out) != void_ptr(&B)) )
+        {
+        arma_debug_print("glue_times_redirect<2>::apply(): inv tiny matrix optimisation");
+        
+        Mat<eT> AA(N, N, arma_nozeros_indicator());
+        
+        arrayops::copy(AA.memptr(), A.memptr(), AA.n_elem);
+        
+        bool inv_status = false;
+        
+        if(N == 1)  { const eT a = AA[0]; AA[0] = eT(1) / a; inv_status = (a != eT(0)); }
+        if(N == 2)  { inv_status = op_inv_gen_full::apply_tiny_2x2(AA); }
+        if(N == 3)  { inv_status = op_inv_gen_full::apply_tiny_3x3(AA); }
+        
+        if(inv_status)  { glue_times::apply<eT,false,false,false>(out, AA, B, eT(0)); return; }
+        
+        arma_debug_print("glue_times_redirect<2>::apply(): inv tiny matrix optimisation failed");
+        
+        // fallthrough if optimisation failed
+        }
+      }
+    
    Mat<eT> A = A_strip.M;
    
    arma_conform_check( (A.is_square() == false), "inv(): given matrix must be square sized" );
@@ -2361,7 +2361,7 @@ gmm_diag<eT>::em_iterate(const Mat<eT>& X, const uword max_iter, const eT var_fl
      get_cout_stream().flush();
      }
    
-    if(arma_isfinite(new_avg_log_p) == false)  { return false; }
+    if(arma_isnonfinite(new_avg_log_p))  { return false; }
    
    if(std::abs(old_avg_log_p - new_avg_log_p) <= Datum<eT>::eps)  { break; }
    
@@ -2473,7 +2473,7 @@ gmm_diag<eT>::em_update_params
    {
    const eT acc_norm_lhood = (std::max)( final_acc_norm_lhoods[g], std::numeric_limits<eT>::min() );
    
-    if(arma_isfinite(acc_norm_lhood) == false)  { continue; }
+    if(arma_isnonfinite(acc_norm_lhood))  { continue; }
    
    eT* acc_mean_mem = final_acc_means.colptr(g);
    eT* acc_dcov_mem = final_acc_dcovs.colptr(g);
@@ -2488,7 +2488,7 @@ gmm_diag<eT>::em_update_params
      acc_mean_mem[d] = tmp1;
      acc_dcov_mem[d] = tmp2;
      
-      if(arma_isfinite(tmp2) == false)  { ok = false; }
+      if(arma_isnonfinite(tmp2))  { ok = false; }
      }
    
    
@@ -2392,7 +2392,7 @@ gmm_full<eT>::em_iterate(const Mat<eT>& X, const uword max_iter, const eT var_fl
      get_cout_stream().flush();
      }
    
-    if(arma_isfinite(new_avg_log_p) == false)  { return false; }
+    if(arma_isnonfinite(new_avg_log_p))  { return false; }
    
    if(std::abs(old_avg_log_p - new_avg_log_p) <= Datum<eT>::eps)  { break; }
    
@@ -2516,7 +2516,7 @@ gmm_full<eT>::em_update_params
    {
    const eT acc_norm_lhood = (std::max)( final_acc_norm_lhoods[g], std::numeric_limits<eT>::min() );
    
-    if(arma_isfinite(acc_norm_lhood) == false)  { continue; }
+    if(arma_isnonfinite(acc_norm_lhood))  { continue; }
    
    eT* acc_mean_mem = final_acc_means.colptr(g);
    
@@ -16,7 +16,7 @@
 // 
 // ------------------------------------------------------------------------
 // 
-// This file includes portions of SuperLU 5.2 software,
+// This file includes portions of SuperLU 7.0 software,
 // licensed under the following conditions.
 // 
 // Copyright (c) 2003, The Regents of the University of California, through
@@ -64,16 +64,23 @@
 // and manually specify a few SuperLU structures and function prototypes.
 //
 // CAVEAT:
-// This code requires SuperLU version 5.2,
-// and assumes that newer 5.x versions will have no API changes.
+// This code requires SuperLU version 7.0, and assumes that newer 7.x versions have no API changes.

 namespace arma
 {
 namespace superlu
  {
-  // slu_*defs.h has int typedefed to int_t.
-  // I'll just write it as int for simplicity, where I can, but supermatrix.h needs int_t.
-  typedef int int_t;
+  // superlu_config.h uses either int or int64_t as int_t
+  
+  #if defined(ARMA_SUPERLU_64BIT_INT)
+    #if defined(INT64_MAX)
+      typedef std::int64_t int_t;
+    #else
+      typedef long long    int_t;
+    #endif
+  #else
+      typedef int          int_t;
+  #endif
  }
 }

@@ -84,7 +91,7 @@ namespace arma
 namespace superlu
  {
  // Include supermatrix.h.  This gives us SuperMatrix.
-  // Put it in the slu namespace.
+  // Put it in the superlu namespace.
  // For versions of SuperLU I am familiar with, supermatrix.h does not include any other files.
  // Therefore, putting it in the superlu namespace is reasonably safe.
  // This same reasoning is true for superlu_enum_consts.h.
@@ -120,7 +127,7 @@ namespace superlu
  
  #undef ARMA_SLU_STR1
  #undef ARMA_SLU_STR2
-    
+  
  #undef ARMA_SLU_HEADER_A
  #undef ARMA_SLU_HEADER_B
  
@@ -130,7 +137,7 @@ namespace superlu
      {
      int*    panel_histo;
      double* utime;
-      float*  ops;
+      float*  ops;        // NOTE: orig definition is flops_t* ops, where flops_t = float
      int     TinyPivots;
      int     RefineSteps;
      int     expansions;
@@ -174,16 +181,16 @@ namespace superlu
    
    typedef struct e_node
      {
-      int   size;
+      int_t size;
      void* mem;
      } ExpHeader;
    
    typedef struct
      {
-      int   size;
-      int   used;
-      int   top1;
-      int   top2;
+      int_t size;
+      int_t used;
+      int_t top1;
+      int_t top2;
      void* array;
      } LU_stack_t;
    
@@ -191,16 +198,16 @@ namespace superlu
      {
      int*       xsup;
      int*       supno;   
-      int*       lsub;
-      int*       xlsub;
+      int_t*     lsub;
+      int_t*     xlsub;
      void*      lusup;
-      int*       xlusup;
+      int_t*     xlusup;
      void*      ucol;
-      int*       usub;
-      int*       xusub;
-      int        nzlmax;
-      int        nzumax;
-      int        nzlumax;
+      int_t*     usub;
+      int_t*     xusub;
+      int_t      nzlmax;
+      int_t      nzumax;
+      int_t      nzlumax;
      int        n;
      LU_space_t MemModel;
      int        num_expansions;
@@ -283,23 +290,23 @@ namespace superlu
    {
    int*    panel_histo;
    double* utime;
-    float*  ops;
+    float*  ops;        // NOTE: orig definition is flops_t* ops, where flops_t = float
    int     TinyPivots;
    int     RefineSteps;
    int     expansions;
    } SuperLUStat_t;
  
-  typedef enum {NO, YES}                                          yes_no_t;
+  typedef enum {NO, YES}                                                yes_no_t;
  typedef enum {DOFACT, SamePattern, SamePattern_SameRowPerm, FACTORED} fact_t;
-  typedef enum {NOROWPERM, LargeDiag, MY_PERMR}                   rowperm_t;
+  typedef enum {NOROWPERM, LargeDiag_MC64, LargeDiag_HWPM, MY_PERMR}    rowperm_t;
  typedef enum {NATURAL, MMD_ATA, MMD_AT_PLUS_A, COLAMD,
-                METIS_AT_PLUS_A, PARMETIS, ZOLTAN, MY_PERMC}      colperm_t;
-  typedef enum {NOTRANS, TRANS, CONJ}                             trans_t;
-  typedef enum {NOREFINE, SLU_SINGLE=1, SLU_DOUBLE, SLU_EXTRA}    IterRefine_t;
-  typedef enum {SYSTEM, USER}                                     LU_space_t;
-  typedef enum {ONE_NORM, TWO_NORM, INF_NORM}                     norm_t;
-  typedef enum {SILU, SMILU_1, SMILU_2, SMILU_3}                  milu_t;
-
+                METIS_AT_PLUS_A, PARMETIS, METIS_ATA, ZOLTAN, MY_PERMC} colperm_t;
+  typedef enum {NOTRANS, TRANS, CONJ}                                   trans_t;
+  typedef enum {NOREFINE, SLU_SINGLE=1, SLU_DOUBLE, SLU_EXTRA}          IterRefine_t;
+  typedef enum {SYSTEM, USER}                                           LU_space_t;
+  typedef enum {ONE_NORM, TWO_NORM, INF_NORM}			        norm_t;
+  typedef enum {SILU, SMILU_1, SMILU_2, SMILU_3}			milu_t;
+  
  typedef struct
    {
    fact_t        Fact;
@@ -352,16 +359,16 @@ namespace superlu
  
  typedef struct e_node
    {
-    int   size;
+    int_t size;
    void* mem;
    } ExpHeader;
  
  typedef struct
    {
-    int   size;
-    int   used;
-    int   top1;
-    int   top2;
+    int_t size;
+    int_t used;
+    int_t top1;
+    int_t top2;
    void* array;
    } LU_stack_t;
  
@@ -369,16 +376,16 @@ namespace superlu
    {
    int*       xsup;
    int*       supno;   
-    int*       lsub;
-    int*       xlsub;
+    int_t*     lsub;
+    int_t*     xlsub;
    void*      lusup;
-    int*       xlusup;
+    int_t*     xlusup;
    void*      ucol;
-    int*       usub;
-    int*       xusub;
-    int        nzlmax;
-    int        nzumax;
-    int        nzlumax;
+    int_t*     usub;
+    int_t*     xusub;
+    int_t      nzlmax;
+    int_t      nzumax;
+    int_t      nzlumax;
    int        n;
    LU_space_t MemModel;
    int        num_expansions;
@@ -21,6 +21,9 @@



+struct mtSpOp_dual_aux_indicator {};
+
+
 template<typename out_eT, typename T1, typename spop_type>
 class mtSpOp : public SpBase< out_eT, mtSpOp<out_eT, T1, spop_type> >
  {
@@ -39,6 +42,7 @@ class mtSpOp : public SpBase< out_eT, mtSpOp<out_eT, T1, spop_type> >
  inline           mtSpOp(const T1& in_m, const in_eT in_aux);
  inline           mtSpOp(const T1& in_m, const uword aux_uword_a, const uword aux_uword_b);
  inline           mtSpOp(const char junk, const T1& in_m, const out_eT in_aux);
+  inline           mtSpOp(const mtSpOp_dual_aux_indicator&, const T1& in_m, const in_eT in_aux_a, const out_eT in_aux_b);
  inline          ~mtSpOp();
  
  template<typename eT2>
@@ -67,6 +67,18 @@ mtSpOp<out_eT, T1, op_type>::mtSpOp(const char junk, const T1& in_m, const out_e



+template<typename out_eT, typename T1, typename op_type>
+inline
+mtSpOp<out_eT, T1, op_type>::mtSpOp(const mtSpOp_dual_aux_indicator&, const T1& in_m, const typename T1::elem_type in_aux_a, const out_eT in_aux_b)
+  : m         (in_m    )
+  , aux       (in_aux_a)
+  , aux_out_eT(in_aux_b)
+  {
+  arma_debug_sigprint();
+  }
+
+
+
 template<typename out_eT, typename T1, typename op_type>
 inline
 mtSpOp<out_eT, T1, op_type>::~mtSpOp()
@@ -104,15 +104,29 @@ op_expmat::apply_direct(Mat<typename T1::elem_type>& out, const Base<typename T1
    return true;
    }
  
+  // trace reduction
+  
+  const eT     diag_shift = arma::trace(A) / T(A.n_rows);
+  const eT exp_diag_shift = std::exp(diag_shift);
+  
+  const bool do_trace_reduction = arma_isfinite(diag_shift) && arma_isfinite(exp_diag_shift) && (exp_diag_shift != eT(0)) && ( (is_cx<eT>::yes) ? (std::abs(diag_shift) > T(0)) : (access::tmp_real(diag_shift) > T(0)) );
+  
+  if(do_trace_reduction)
+    {
+    arma_debug_print("op_expmat: diag_shift: ", diag_shift);
+    
+    A.diag() -= diag_shift;
+    }
+  
  const T norm_val = arma::norm(A, "inf");
  
-  if(arma_isfinite(norm_val) == false)  { return false; }
+  if(arma_isnonfinite(norm_val))  { return false; }
  
-  const double log2_val = (norm_val > T(0)) ? double(eop_aux::log2(norm_val)) : double(0);
+  int exponent = int(0);  std::frexp(norm_val, &exponent);
  
-  int exponent = int(0);  std::frexp(log2_val, &exponent);
+  const uword s = (std::min)( uword( (std::max)(int(0), exponent) ), uword(1023) );
  
-  const uword s = uword( (std::max)(int(0), exponent + int(1)) );
+  arma_debug_print("op_expmat: s: ", s);
  
  A /= eT(eop_aux::pow(double(2), double(s)));
  
@@ -125,7 +139,7 @@ op_expmat::apply_direct(Mat<typename T1::elem_type>& out, const Base<typename T1
  
  bool positive = true;
  
-  const uword N = 6;
+  const uword N = 8;
  
  for(uword i = 2; i <= N; ++i)
    {
@@ -148,6 +162,9 @@ op_expmat::apply_direct(Mat<typename T1::elem_type>& out, const Base<typename T1
  
  for(uword i=0; i < s; ++i)  { out = out * out; }
  
+  // inverse trace reduction
+  if(do_trace_reduction)  { out *= exp_diag_shift; }
+  
  return true;
  }

@@ -127,4 +127,15 @@ class op_find_nan



+class op_find_nonnan
+  : public traits_op_col
+  {
+  public:
+
+  template<typename T1>
+  inline static void apply(Mat<uword>& out, const mtOp<uword, T1, op_find_nonnan>& X);
+  };
+
+
+
 //! @}
@@ -585,7 +585,7 @@ op_find_nonfinite::apply(Mat<uword>& out, const mtOp<uword, T1, op_find_nonfinit
    
    for(uword i=0; i<n_elem; ++i)
      {
-      if( arma_isfinite(Pea[i]) == false )  { indices_mem[count] = i; count++; }
+      if( arma_isnonfinite(Pea[i]) )  { indices_mem[count] = i; count++; }
      }
    }
  else
@@ -598,7 +598,7 @@ op_find_nonfinite::apply(Mat<uword>& out, const mtOp<uword, T1, op_find_nonfinit
    for(uword col=0; col<n_cols; ++col)
    for(uword row=0; row<n_rows; ++row)
      {
-      if( arma_isfinite(P.at(row,col)) == false )  { indices_mem[count] = i; count++; }
+      if( arma_isnonfinite(P.at(row,col)) )  { indices_mem[count] = i; count++; }
      
      i++;
      }
@@ -657,4 +657,52 @@ op_find_nan::apply(Mat<uword>& out, const mtOp<uword, T1, op_find_nan>& X)



+template<typename T1>
+inline
+void
+op_find_nonnan::apply(Mat<uword>& out, const mtOp<uword, T1, op_find_nonnan>& X)
+  {
+  arma_debug_sigprint();
+  
+  if(arma_config::fast_math_warn)  { arma_warn(1, "find_nonnan(): detection of non-finite values is not reliable in fast math mode"); }
+  
+  const Proxy<T1> P(X.m);
+  
+  const uword n_elem = P.get_n_elem();
+  
+  Mat<uword> indices(n_elem, 1, arma_nozeros_indicator());
+  
+  uword* indices_mem = indices.memptr();
+  uword  count       = 0;
+  
+  if(Proxy<T1>::use_at == false)
+    {
+    const typename Proxy<T1>::ea_type Pea = P.get_ea();
+    
+    for(uword i=0; i < n_elem; ++i)
+      {
+      if( arma_isnan(Pea[i]) == false )  { indices_mem[count] = i; ++count; }
+      }
+    }
+  else
+    {
+    const uword n_rows = P.get_n_rows();
+    const uword n_cols = P.get_n_cols();
+    
+    uword i = 0;
+    
+    for(uword col=0; col < n_cols; ++col)
+    for(uword row=0; row < n_rows; ++row)
+      {
+      if( arma_isnan(P.at(row,col)) == false )  { indices_mem[count] = i; ++count; }
+      
+      ++i;
+      }
+    }
+  
+  out.steal_mem_col(indices, count);
+  }
+
+
+
 //! @}
@@ -66,8 +66,8 @@ op_hist::apply_noalias(Mat<uword>& out, const Mat<eT>& A, const uword n_bins, co
    max_val += (n_bins/2);
    }
  
-  if(arma_isfinite(min_val) == false) { min_val = priv::most_neg<eT>(); }
-  if(arma_isfinite(max_val) == false) { max_val = priv::most_pos<eT>(); }
+  if(arma_isnonfinite(min_val)) { min_val = priv::most_neg<eT>(); }
+  if(arma_isnonfinite(max_val)) { max_val = priv::most_pos<eT>(); }
  
  Col<eT> c(n_bins, arma_nozeros_indicator());
  eT* c_mem = c.memptr();
@@ -31,14 +31,8 @@ class 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_noalias(Mat<typename T1::elem_type>& out, const Proxy<T1>& P, const uword dim);
-  
-  template<typename T1>
-  inline static void apply_noalias_unwrap(Mat<typename T1::elem_type>& out, const Proxy<T1>& P, const uword dim);
-  
-  template<typename T1>
-  inline static void apply_noalias_proxy(Mat<typename T1::elem_type>& out, const Proxy<T1>& P, const uword dim);
+  template<typename eT>
+  inline static void apply_noalias(Mat<eT>& out, const Mat<eT>& X, const uword dim);
  
  
  // cubes
@@ -46,60 +40,28 @@ class op_mean
  template<typename T1>
  inline static void apply(Cube<typename T1::elem_type>& out, const OpCube<T1,op_mean>& in);
  
-  template<typename T1>
-  inline static void apply_noalias(Cube<typename T1::elem_type>& out, const ProxyCube<T1>& P, const uword dim);
-  
-  template<typename T1>
-  inline static void apply_noalias_unwrap(Cube<typename T1::elem_type>& out, const ProxyCube<T1>& P, const uword dim);
-  
-  template<typename T1>
-  inline static void apply_noalias_proxy(Cube<typename T1::elem_type>& out, const ProxyCube<T1>& P, const uword dim);
-  
+  template<typename eT>
+  inline static void apply_noalias(Cube<eT>& out, const Cube<eT>& X, const uword dim);
  
  //
  
  template<typename eT>
-  inline static eT direct_mean(const eT* const X, const uword N);
+  inline static eT direct_mean(const eT* X_mem, const uword N);
  
  template<typename eT>
-  inline static eT direct_mean_robust(const eT* const X, const uword N);
-  
-
-  //
-  
-  template<typename eT>
-  inline static eT direct_mean(const Mat<eT>& X, const uword row);
-  
-  template<typename eT>
-  inline static eT direct_mean_robust(const Mat<eT>& X, const uword row);
-  
-
-  //
-  
-  template<typename eT>
-  inline static eT mean_all(const subview<eT>& X);
-  
-  template<typename eT>
-  inline static eT mean_all_robust(const subview<eT>& X);
-  
-  
-  //
-  
-  template<typename eT>
-  inline static eT mean_all(const diagview<eT>& X);
-  
-  template<typename eT>
-  inline static eT mean_all_robust(const diagview<eT>& X);
+  inline static eT direct_mean_robust(const eT old_mean, const eT* X_mem, const uword N);
  
  
  //
  
  template<typename T1>
-  inline static typename T1::elem_type mean_all(const Op<T1,op_vectorise_col>& X);
+  inline static typename T1::elem_type mean_all(const T1& X);
  
  template<typename T1>
-  inline static typename T1::elem_type mean_all(const Base<typename T1::elem_type, T1>& X);
+  inline static typename T1::elem_type mean_all(const Op<T1, op_omit>& X);
  
+  template<typename eT, typename functor>
+  inline static eT mean_all_omit(const eT* X_mem, const uword N, functor is_omitted);
  
  //
  
@@ -111,5 +73,4 @@ class op_mean
  };


-
 //! @}
@@ -31,60 +31,35 @@ op_mean::apply(Mat<typename T1::elem_type>& out, const Op<T1,op_mean>& in)
  typedef typename T1::elem_type eT;
  
  const uword dim = in.aux_uword_a;
+  
  arma_conform_check( (dim > 1), "mean(): parameter 'dim' must be 0 or 1" );
  
-  const Proxy<T1> P(in.m);
+  const quasi_unwrap<T1> U(in.m);
  
-  if(P.is_alias(out) == false)
-    {
-    op_mean::apply_noalias(out, P, dim);
-    }
-  else
+  if(U.is_alias(out))
    {
    Mat<eT> tmp;
    
-    op_mean::apply_noalias(tmp, P, dim);
+    op_mean::apply_noalias(tmp, U.M, dim);
    
    out.steal_mem(tmp);
    }
-  }
-
-
-
-template<typename T1>
-inline
-void
-op_mean::apply_noalias(Mat<typename T1::elem_type>& out, const Proxy<T1>& P, const uword dim)
-  {
-  arma_debug_sigprint();
-  
-  if((is_Mat<typename Proxy<T1>::stored_type>::value) || (arma_config::openmp && Proxy<T1>::use_mp))
-    {
-    op_mean::apply_noalias_unwrap(out, P, dim);
-    }
  else
    {
-    op_mean::apply_noalias_proxy(out, P, dim);
+    op_mean::apply_noalias(out, U.M, dim);
    }
  }



-template<typename T1>
+template<typename eT>
 inline
 void
-op_mean::apply_noalias_unwrap(Mat<typename T1::elem_type>& out, const Proxy<T1>& P, const uword dim)
+op_mean::apply_noalias(Mat<eT>& out, const Mat<eT>& X, const uword dim)
  {
  arma_debug_sigprint();
  
-  typedef typename T1::elem_type            eT;
-  typedef typename get_pod_type<eT>::result  T;
-  
-  typedef typename Proxy<T1>::stored_type P_stored_type;
-  
-  const unwrap<P_stored_type> tmp(P.Q);
-  
-  const typename unwrap<P_stored_type>::stored_type& X = tmp.M;
+  typedef typename get_pod_type<eT>::result T;
  
  const uword X_n_rows = X.n_rows;
  const uword X_n_cols = X.n_cols;
@@ -113,21 +88,25 @@ op_mean::apply_noalias_unwrap(Mat<typename T1::elem_type>& out, const Proxy<T1>&
    
    for(uword col=0; col < X_n_cols; ++col)
      {
-      const eT* col_mem = X.colptr(col);
-      
-      for(uword row=0; row < X_n_rows; ++row)
-        {
-        out_mem[row] += col_mem[row];
-        }
+      arrayops::inplace_plus(out_mem, X.colptr(col), X_n_rows);
      }
    
    out /= T(X_n_cols);
    
-    for(uword row=0; row < X_n_rows; ++row)
+    if(out.internal_has_nonfinite())
      {
-      if(arma_isfinite(out_mem[row]) == false)
+      podarray<eT> tmp;
+      
+      for(uword row=0; row < X_n_rows; ++row)
        {
-        out_mem[row] = op_mean::direct_mean_robust( X, row );
+        const eT old_mean = out_mem[row];
+        
+        if(arma_isnonfinite(old_mean))
+          {
+          tmp.copy_row(X, row);
+          
+          out_mem[row] = op_mean::direct_mean_robust(old_mean, tmp.memptr(), tmp.n_elem);
+          }
        }
      }
    }
@@ -135,76 +114,7 @@ op_mean::apply_noalias_unwrap(Mat<typename T1::elem_type>& out, const Proxy<T1>&



-template<typename T1>
-inline
-void
-op_mean::apply_noalias_proxy(Mat<typename T1::elem_type>& out, const Proxy<T1>& P, const uword dim)
-  {
-  arma_debug_sigprint();
-  
-  typedef typename T1::elem_type            eT;
-  typedef typename get_pod_type<eT>::result  T;
-  
-  const uword P_n_rows = P.get_n_rows();
-  const uword P_n_cols = P.get_n_cols();
-  
-  if(dim == 0)
-    {
-    out.set_size((P_n_rows > 0) ? 1 : 0, P_n_cols);
-    
-    if(P_n_rows == 0)  { return; }
-    
-    eT* out_mem = out.memptr();
-    
-    for(uword col=0; col < P_n_cols; ++col)
-      {
-      eT val1 = eT(0);
-      eT val2 = eT(0);
-      
-      uword i,j;
-      for(i=0, j=1; j < P_n_rows; i+=2, j+=2)
-        {
-        val1 += P.at(i,col);
-        val2 += P.at(j,col);
-        }
-      
-      if(i < P_n_rows)
-        {
-        val1 += P.at(i,col);
-        }
-      
-      out_mem[col] = (val1 + val2) / T(P_n_rows);
-      }
-    }
-  else
-  if(dim == 1)
-    {
-    out.zeros(P_n_rows, (P_n_cols > 0) ? 1 : 0);
-    
-    if(P_n_cols == 0)  { return; }
-    
-    eT* out_mem = out.memptr();
-    
-    for(uword col=0; col < P_n_cols; ++col)
-    for(uword row=0; row < P_n_rows; ++row)
-      {
-      out_mem[row] += P.at(row,col);
-      }
-    
-    out /= T(P_n_cols);
-    }
-  
-  if(out.internal_has_nonfinite())
-    {
-    // TODO: replace with dedicated handling to avoid unwrapping
-    op_mean::apply_noalias_unwrap(out, P, dim);
-    }
-  }
-
-
-
 //
-// cubes



@@ -218,60 +128,35 @@ op_mean::apply(Cube<typename T1::elem_type>& out, const OpCube<T1,op_mean>& in)
  typedef typename T1::elem_type eT;
  
  const uword dim = in.aux_uword_a;
+  
  arma_conform_check( (dim > 2), "mean(): parameter 'dim' must be 0 or 1 or 2" );
  
-  const ProxyCube<T1> P(in.m);
+  const unwrap_cube<T1> U(in.m);
  
-  if(P.is_alias(out) == false)
-    {
-    op_mean::apply_noalias(out, P, dim);
-    }
-  else
+  if(U.is_alias(out))
    {
    Cube<eT> tmp;
    
-    op_mean::apply_noalias(tmp, P, dim);
+    op_mean::apply_noalias(tmp, U.M, dim);
    
    out.steal_mem(tmp);
    }
-  }
-
-
-
-template<typename T1>
-inline
-void
-op_mean::apply_noalias(Cube<typename T1::elem_type>& out, const ProxyCube<T1>& P, const uword dim)
-  {
-  arma_debug_sigprint();
-  
-  if((is_Cube<typename ProxyCube<T1>::stored_type>::value) || (arma_config::openmp && ProxyCube<T1>::use_mp))
-    {
-    op_mean::apply_noalias_unwrap(out, P, dim);
-    }
  else
    {
-    op_mean::apply_noalias_proxy(out, P, dim);
+    op_mean::apply_noalias(out, U.M, dim);
    }
  }



-template<typename T1>
+template<typename eT>
 inline
 void
-op_mean::apply_noalias_unwrap(Cube<typename T1::elem_type>& out, const ProxyCube<T1>& P, const uword dim)
+op_mean::apply_noalias(Cube<eT>& out, const Cube<eT>& X, const uword dim)
  {
  arma_debug_sigprint();
  
-  typedef typename T1::elem_type            eT;
-  typedef typename get_pod_type<eT>::result  T;
-  
-  typedef typename ProxyCube<T1>::stored_type P_stored_type;
-  
-  const unwrap_cube<P_stored_type> U(P.Q);
-  
-  const Cube<eT>& X = U.M;
+  typedef typename get_pod_type<eT>::result T;
  
  const uword X_n_rows   = X.n_rows;
  const uword X_n_cols   = X.n_cols;
@@ -306,23 +191,27 @@ op_mean::apply_noalias_unwrap(Cube<typename T1::elem_type>& out, const ProxyCube
      
      for(uword col=0; col < X_n_cols; ++col)
        {
-        const eT* col_mem = X.slice_colptr(slice,col);
+        arrayops::inplace_plus(out_mem, X.slice_colptr(slice,col), X_n_rows);
+        }
+      
+      for(uword row=0; row < X_n_rows; ++row)  { out_mem[row] /= T(X_n_cols); }
+      
+      if(arrayops::is_finite(out_mem, X_n_rows) == false)
+        {
+        const Mat<eT> tmp_mat('j', X.slice_memptr(slice), X_n_rows, X_n_cols);
+        
+        podarray<eT> tmp_vec;
        
        for(uword row=0; row < X_n_rows; ++row)
          {
-          out_mem[row] += col_mem[row];
-          }
-        }
-      
-      const Mat<eT> tmp('j', X.slice_memptr(slice), X_n_rows, X_n_cols);
-      
-      for(uword row=0; row < X_n_rows; ++row)
-        {
-        out_mem[row] /= T(X_n_cols);
-        
-        if(arma_isfinite(out_mem[row]) == false)
-          {
-          out_mem[row] = op_mean::direct_mean_robust( tmp, row );
+          const eT old_mean = out_mem[row];
+          
+          if(arma_isnonfinite(old_mean))
+            {
+            tmp_vec.copy_row(tmp_mat, row);
+            
+            out_mem[row] = op_mean::direct_mean_robust(old_mean, tmp_vec.memptr(), tmp_vec.n_elem);
+            }
          }
        }
      }
@@ -343,19 +232,21 @@ op_mean::apply_noalias_unwrap(Cube<typename T1::elem_type>& out, const ProxyCube
    
    out /= T(X_n_slices);
    
-    podarray<eT> tmp(X_n_slices);
-      
-    for(uword col=0; col < X_n_cols; ++col)
-    for(uword row=0; row < X_n_rows; ++row)
+    if(out.internal_has_nonfinite())
      {
-      if(arma_isfinite(out.at(row,col,0)) == false)
+      podarray<eT> tmp(X_n_slices);
+      
+      for(uword col=0; col < X_n_cols; ++col)
+      for(uword row=0; row < X_n_rows; ++row)
        {
-        for(uword slice=0; slice < X_n_slices; ++slice)
-          {
-          tmp[slice] = X.at(row,col,slice);
-          }
+        const eT old_mean = out.at(row,col,0);
        
-        out.at(row,col,0) = op_mean::direct_mean_robust(tmp.memptr(), X_n_slices);
+        if(arma_isnonfinite(old_mean))
+          {
+          for(uword slice=0; slice < X_n_slices; ++slice)  { tmp[slice] = X.at(row,col,slice); }
+          
+          out.at(row,col,0) = op_mean::direct_mean_robust(old_mean, tmp.memptr(), tmp.n_elem);
+          }
        }
      }
    }
@@ -363,21 +254,6 @@ op_mean::apply_noalias_unwrap(Cube<typename T1::elem_type>& out, const ProxyCube



-template<typename T1>
-inline
-void
-op_mean::apply_noalias_proxy(Cube<typename T1::elem_type>& out, const ProxyCube<T1>& P, const uword dim)
-  {
-  arma_debug_sigprint();
-  
-  op_mean::apply_noalias_unwrap(out, P, dim);
-  
-  // TODO: implement specialised handling
-  }
-
-
-
-
 //


@@ -385,15 +261,15 @@ op_mean::apply_noalias_proxy(Cube<typename T1::elem_type>& out, const ProxyCube<
 template<typename eT>
 inline
 eT
-op_mean::direct_mean(const eT* const X, const uword n_elem)
+op_mean::direct_mean(const eT* X_mem, const uword N)
  {
  arma_debug_sigprint();
  
  typedef typename get_pod_type<eT>::result T;
  
-  const eT result = arrayops::accumulate(X, n_elem) / T(n_elem);
+  const eT mean = arrayops::accumulate(X_mem, N) / T(N);
  
-  return arma_isfinite(result) ? result : op_mean::direct_mean_robust(X, n_elem);
+  return arma_isfinite(mean) ? mean : op_mean::direct_mean_robust(mean, X_mem, N);
  }


@@ -401,7 +277,7 @@ op_mean::direct_mean(const eT* const X, const uword n_elem)
 template<typename eT>
 inline
 eT
-op_mean::direct_mean_robust(const eT* const X, const uword n_elem)
+op_mean::direct_mean_robust(const eT old_mean, const eT* X_mem, const uword N)
  {
  arma_debug_sigprint();
  
@@ -409,25 +285,13 @@ op_mean::direct_mean_robust(const eT* const X, const uword n_elem)
  
  typedef typename get_pod_type<eT>::result T;
  
-  uword i,j;
+  if(arrayops::is_finite(X_mem, N) == false)  { return old_mean; }
  
  eT r_mean = eT(0);
  
-  for(i=0, j=1; j<n_elem; i+=2, j+=2)
+  for(uword i=0; i < N; ++i)
    {
-    const eT Xi = X[i];
-    const eT Xj = X[j];
-    
-    r_mean = r_mean + (Xi - r_mean)/T(j);    // we need i+1, and j is equivalent to i+1 here
-    r_mean = r_mean + (Xj - r_mean)/T(j+1);
-    }
-  
-  
-  if(i < n_elem)
-    {
-    const eT Xi = X[i];
-    
-    r_mean = r_mean + (Xi - r_mean)/T(i+1);
+    r_mean = r_mean + (X_mem[i] - r_mean) / T(i+1);
    }
  
  return r_mean;
@@ -435,266 +299,128 @@ op_mean::direct_mean_robust(const eT* const X, const uword n_elem)



-template<typename eT>
-inline
-eT
-op_mean::direct_mean(const Mat<eT>& X, const uword row)
-  {
-  arma_debug_sigprint();
-  
-  typedef typename get_pod_type<eT>::result T;
-  
-  const uword X_n_cols = X.n_cols;
-  
-  eT val = eT(0);
-  
-  uword i,j;
-  for(i=0, j=1; j < X_n_cols; i+=2, j+=2)
-    {
-    val += X.at(row,i);
-    val += X.at(row,j);
-    }
-  
-  if(i < X_n_cols)
-    {
-    val += X.at(row,i);
-    }
-  
-  const eT result = val / T(X_n_cols);
-  
-  return arma_isfinite(result) ? result : op_mean::direct_mean_robust(X, row);
-  }
-
-
-
-template<typename eT>
-inline
-eT
-op_mean::direct_mean_robust(const Mat<eT>& X, const uword row)
-  {
-  arma_debug_sigprint();
-  
-  typedef typename get_pod_type<eT>::result T;
-  
-  const uword X_n_cols = X.n_cols;
-  
-  eT r_mean = eT(0);
-  
-  for(uword col=0; col < X_n_cols; ++col)
-    {
-    r_mean = r_mean + (X.at(row,col) - r_mean)/T(col+1);
-    }
-  
-  return r_mean;
-  }
-
-
-
-template<typename eT>
-inline
-eT
-op_mean::mean_all(const subview<eT>& X)
-  {
-  arma_debug_sigprint();
-  
-  typedef typename get_pod_type<eT>::result T;
-  
-  const uword X_n_rows = X.n_rows;
-  const uword X_n_cols = X.n_cols;
-  const uword X_n_elem = X.n_elem;
-  
-  if(X_n_elem == 0)
-    {
-    arma_conform_check(true, "mean(): object has no elements");
-    
-    return Datum<eT>::nan;
-    }
-  
-  eT val = eT(0);
-  
-  if(X_n_rows == 1)
-    {
-    const Mat<eT>& A = X.m;
-    
-    const uword start_row = X.aux_row1;
-    const uword start_col = X.aux_col1;
-    
-    const uword end_col_p1 = start_col + X_n_cols;
-    
-    uword i,j;
-    for(i=start_col, j=start_col+1; j < end_col_p1; i+=2, j+=2)
-      {
-      val += A.at(start_row, i);
-      val += A.at(start_row, j);
-      }
-    
-    if(i < end_col_p1)
-      {
-      val += A.at(start_row, i);
-      }
-    }
-  else
-    {
-    for(uword col=0; col < X_n_cols; ++col)
-      {
-      val += arrayops::accumulate(X.colptr(col), X_n_rows);
-      }
-    }
-  
-  const eT result = val / T(X_n_elem);
-  
-  return arma_isfinite(result) ? result : op_mean::mean_all_robust(X);
-  }
-
-
-
-template<typename eT>
-inline 
-eT
-op_mean::mean_all_robust(const subview<eT>& X)
-  {
-  arma_debug_sigprint();
-  
-  typedef typename get_pod_type<eT>::result T;
-  
-  const uword X_n_rows = X.n_rows;
-  const uword X_n_cols = X.n_cols;
-  
-  const uword start_row = X.aux_row1;
-  const uword start_col = X.aux_col1;
-  
-  const uword end_row_p1 = start_row + X_n_rows;
-  const uword end_col_p1 = start_col + X_n_cols;
-  
-  const Mat<eT>& A = X.m;
-  
-  
-  eT r_mean = eT(0);
-  
-  if(X_n_rows == 1)
-    {
-    uword i=0;
-    
-    for(uword col = start_col; col < end_col_p1; ++col, ++i)
-      {
-      r_mean = r_mean + (A.at(start_row,col) - r_mean)/T(i+1);
-      }
-    }
-  else
-    {
-    uword i=0;
-    
-    for(uword col = start_col; col < end_col_p1; ++col)
-    for(uword row = start_row; row < end_row_p1; ++row, ++i)
-      {
-      r_mean = r_mean + (A.at(row,col) - r_mean)/T(i+1);
-      }
-    }
-  
-  return r_mean;
-  }
-
-
-
-template<typename eT>
-inline 
-eT
-op_mean::mean_all(const diagview<eT>& X)
-  {
-  arma_debug_sigprint();
-  
-  typedef typename get_pod_type<eT>::result T;
-  
-  const uword X_n_elem = X.n_elem;
-  
-  if(X_n_elem == 0)
-    {
-    arma_conform_check(true, "mean(): object has no elements");
-    
-    return Datum<eT>::nan;
-    }
-  
-  eT val = eT(0);
-  
-  for(uword i=0; i<X_n_elem; ++i)
-    {
-    val += X[i];
-    }
-  
-  const eT result = val / T(X_n_elem);
-  
-  return arma_isfinite(result) ? result : op_mean::mean_all_robust(X);
-  }
-
-
-
-template<typename eT>
-inline 
-eT
-op_mean::mean_all_robust(const diagview<eT>& X)
-  {
-  arma_debug_sigprint();
-  
-  typedef typename get_pod_type<eT>::result T;
-  
-  const uword X_n_elem = X.n_elem;
-  
-  eT r_mean = eT(0);
-  
-  for(uword i=0; i<X_n_elem; ++i)
-    {
-    r_mean = r_mean + (X[i] - r_mean)/T(i+1);
-    }
-  
-  return r_mean;
-  }
+//



 template<typename T1>
 inline
-typename T1::elem_type 
-op_mean::mean_all(const Op<T1,op_vectorise_col>& X)
-  {
-  arma_debug_sigprint();
-  
-  return op_mean::mean_all(X.m);
-  }
-
-
-
-template<typename T1>
-inline
-typename T1::elem_type 
-op_mean::mean_all(const Base<typename T1::elem_type, T1>& X)
+typename T1::elem_type
+op_mean::mean_all(const T1& X)
  {
  arma_debug_sigprint();
  
  typedef typename T1::elem_type eT;
  
-  const quasi_unwrap<T1> tmp(X.get_ref());
-  const Mat<eT>&     A = tmp.M;
+  const quasi_unwrap<T1> U(X);
  
-  const uword A_n_elem = A.n_elem;
-  
-  if(A_n_elem == 0)
+  if(U.M.n_elem == 0)
    {
    arma_conform_check(true, "mean(): object has no elements");
    
    return Datum<eT>::nan;
    }
  
-  return op_mean::direct_mean(A.memptr(), A_n_elem);
+  return op_mean::direct_mean(U.M.memptr(), U.M.n_elem);
  }



+template<typename T1>
+inline
+typename T1::elem_type
+op_mean::mean_all(const Op<T1, op_omit>& in)
+  {
+  arma_debug_sigprint();
+  
+  typedef typename T1::elem_type eT;
+  
+  const uword omit_mode = in.aux_uword_a;
+  
+  if(arma_config::fast_math_warn)
+    {
+    if(omit_mode == 1)  { arma_warn(1, "omit_nan(): detection of NaN is not reliable in fast math mode"); }
+    if(omit_mode == 2)  { arma_warn(1, "omit_nonfinite(): detection of non-finite values is not reliable in fast math mode"); }
+    }
+  
+  const quasi_unwrap<T1> U(in.m);
+  
+  if(U.M.n_elem == 0)
+    {
+    arma_conform_check(true, "mean(): object has no elements");
+    
+    return Datum<eT>::nan;
+    }
+  
+  auto is_omitted_1 = [](const eT& x) -> bool { return arma_isnan(x);       };
+  auto is_omitted_2 = [](const eT& x) -> bool { return arma_isnonfinite(x); };
+  
+  eT result = eT(0);
+  
+  if(omit_mode == 1)  { result = op_mean::mean_all_omit(U.M.memptr(), U.M.n_elem, is_omitted_1); }
+  if(omit_mode == 2)  { result = op_mean::mean_all_omit(U.M.memptr(), U.M.n_elem, is_omitted_2); }
+  
+  return result;
+  }
+
+
+
+template<typename eT, typename functor>
+inline
+eT
+op_mean::mean_all_omit(const eT* X_mem, const uword N, functor is_omitted)
+  {
+  arma_debug_sigprint();
+  
+  typedef typename get_pod_type<eT>::result T;
+  
+  uword count = 0;
+  eT    acc   = eT(0);
+  
+  for(uword i=0; i < N; ++i)
+    {
+    const eT val = X_mem[i];
+    
+    if(is_omitted(val) == false)  { acc += val; ++count; }
+    }
+  
+  acc /= T(count);
+  
+  if(arma_isfinite(acc))  { return acc; }
+  
+  // handle possible overflow
+  
+  eT r_mean = eT(0);
+  
+  count = 0;
+  
+  for(uword i=0; i < N; ++i)
+    {
+    const eT val = X_mem[i];
+    
+    if(is_omitted(val) == false)
+      {
+      r_mean = r_mean + (val - r_mean) / T(count+1);  // kept as count+1 to use same algorithm as op_mean::direct_mean_robust()
+      
+      ++count;
+      }
+    }
+  
+  return r_mean;
+  }
+
+
+
+//
+
+
+
 template<typename eT>
 arma_inline
 eT
 op_mean::robust_mean(const eT A, const eT B)
  {
-  return A + (B - A)/eT(2);
+  return (arma_isfinite(A) && arma_isfinite(B)) ? eT( A + (B - A)/eT(2) ) : eT( (A+B)/eT(2) );
  }


@@ -704,10 +430,11 @@ arma_inline
 std::complex<T>
 op_mean::robust_mean(const std::complex<T>& A, const std::complex<T>& B)
  {
-  return A + (B - A)/T(2);
+  typedef typename std::complex<T> eT;
+  
+  return (arma_isfinite(A) && arma_isfinite(B)) ? eT( A + (B - A)/T(2) ) : eT( (A+B)/T(2) );
  }



 //! @}
-
@@ -41,7 +41,7 @@ op_median::apply(Mat<typename T1::elem_type>& out, const Op<T1,op_median>& expr)
    {
    Mat<eT> tmp;
    
-    op_median::apply_noalias(out, U.M, dim);
+    op_median::apply_noalias(tmp, U.M, dim);
    
    out.steal_mem(tmp);
    }
@@ -77,4 +77,16 @@ class op_arg



+class op_replace
+  : public traits_op_passthru
+  {
+  public:
+  
+  template<typename eT, typename T1> inline static void apply(Mat<eT>& out, const mtOp<eT,T1,op_replace>& in);
+  
+  template<typename eT, typename T1> inline static void apply(Cube<eT>& out, const mtOpCube<eT,T1,op_replace>& in);
+  };
+
+
+
 //! @}
@@ -411,4 +411,38 @@ op_arg::apply( Cube<typename T1::pod_type>& out, const mtOpCube<typename T1::pod



+template<typename eT, typename T1>
+inline
+void
+op_replace::apply(Mat<eT>& out, const mtOp<eT,T1,op_replace>& in)
+  {
+  arma_debug_sigprint();
+  
+  const eT old_val = in.aux;
+  const eT new_val = in.aux_out_eT;
+  
+  out = in.m;
+  
+  out.replace(old_val, new_val);
+  }
+
+
+
+template<typename eT, typename T1>
+inline
+void
+op_replace::apply(Cube<eT>& out, const mtOpCube<eT,T1,op_replace>& in)
+  {
+  arma_debug_sigprint();
+  
+  const eT old_val = in.aux;
+  const eT new_val = in.aux_out_eT;
+  
+  out = in.m;
+  
+  out.replace(old_val, new_val);
+  }
+
+
+
 //! @}
@@ -31,6 +31,8 @@ op_nonzeros::apply_noalias(Mat<typename T1::elem_type>& out, const Proxy<T1>& P)
  
  typedef typename T1::elem_type eT;
  
+  constexpr eT eT_zero = eT(0);
+  
  const uword N_max = P.get_n_elem();
  
  Mat<eT> tmp(N_max, 1, arma_nozeros_indicator());
@@ -47,7 +49,7 @@ op_nonzeros::apply_noalias(Mat<typename T1::elem_type>& out, const Proxy<T1>& P)
      {
      const eT val = Pea[i];
      
-      if(val != eT(0))  { tmp_mem[N_nz] = val; ++N_nz; }
+      if(val != eT_zero)  { tmp_mem[N_nz] = val; ++N_nz; }
      }
    }
  else
@@ -60,7 +62,7 @@ op_nonzeros::apply_noalias(Mat<typename T1::elem_type>& out, const Proxy<T1>& P)
      {
      const eT val = P.at(row,col);
      
-      if(val != eT(0))  { tmp_mem[N_nz] = val; ++N_nz; }
+      if(val != eT_zero)  { tmp_mem[N_nz] = val; ++N_nz; }
      }
    }
  
@@ -84,11 +86,11 @@ op_nonzeros::apply(Mat<typename T1::elem_type>& out, const Op<T1, op_nonzeros>&
  
  if(P.is_alias(out))
    {
-    Mat<eT> out2;
+    Mat<eT> tmp;
    
-    op_nonzeros::apply_noalias(out2, P);
+    op_nonzeros::apply_noalias(tmp, P);
    
-    out.steal_mem(out2);
+    out.steal_mem(tmp);
    }
  else
    {
@@ -138,7 +138,7 @@ op_norm2est::norm2est
    
    T x_norm = op_norm::vec_norm_2( Proxy< Col<eT> >(x) );
    
-    if(x_norm == T(0) || (arma_isfinite(x_norm) == false) || (x.internal_has_nonfinite()))
+    if( (x_norm == T(0)) || arma_isnonfinite(x_norm) || x.internal_has_nonfinite() )
      {
      randu_filler.fill(x.memptr(), x.n_elem);
      
@@ -155,7 +155,7 @@ op_norm2est::norm2est
    arma_debug_print(arma_str::format("norm2est(): est_old: %e") % est_old);
    arma_debug_print(arma_str::format("norm2est(): est_cur: %e") % est_cur);
    
-    if(arma_isfinite(est_cur) == false)  { return est_old; }
+    if(arma_isnonfinite(est_cur))  { return est_old; }
    
    if( ((std::abs)(est_cur - est_old)) <= (tol * (std::max)(est_cur,est_old)) )  { break; }
    }
@@ -218,7 +218,7 @@ op_norm2est::norm2est
    
    T x_norm = op_norm::vec_norm_2( Proxy< Mat<eT> >(x) );
    
-    if(x_norm == T(0) || (arma_isfinite(x_norm) == false) || (x.internal_has_nonfinite()))
+    if( (x_norm == T(0)) || arma_isnonfinite(x_norm) || x.internal_has_nonfinite() )
      {
      randu_filler.fill(x.memptr(), x.n_elem);
      
@@ -235,7 +235,7 @@ op_norm2est::norm2est
    arma_debug_print(arma_str::format("norm2est(): est_old: %e") % est_old);
    arma_debug_print(arma_str::format("norm2est(): est_cur: %e") % est_cur);
    
-    if(arma_isfinite(est_cur) == false)  { return est_old; }
+    if(arma_isnonfinite(est_cur))  { return est_old; }
    
    if( ((std::abs)(est_cur - est_old)) <= (tol * (std::max)(est_cur,est_old)) )  { break; }
    }
@@ -0,0 +1,50 @@
+// SPDX-License-Identifier: Apache-2.0
+// 
+// Copyright 2008-2016 Conrad Sanderson (http://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
+// http://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 op_omit
+//! @{
+
+
+class op_omit
+  : public traits_op_col
+  {
+  public:
+  
+  template<typename T1>
+  inline static void apply(Mat<typename T1::elem_type>& out, const Op<T1,op_omit>& in);
+  
+  template<typename T1, typename functor>
+  inline static void apply(Mat<typename T1::elem_type>& out, const T1& X, functor is_omitted);
+  };
+
+
+
+class op_omit_cube
+  : public traits_op_col
+  {
+  public:
+  
+  template<typename T1> inline static void apply(Mat<typename T1::elem_type>& out, const CubeToMatOp<T1, op_omit_cube>& in);
+  
+  template<typename T1, typename functor>
+  inline static void apply(Mat<typename T1::elem_type>& out, const T1& X, functor is_omitted);
+  };
+
+
+
+//! @}
@@ -0,0 +1,230 @@
+// SPDX-License-Identifier: Apache-2.0
+// 
+// Copyright 2008-2016 Conrad Sanderson (http://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
+// http://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 op_omit
+//! @{
+
+
+
+template<typename T1>
+inline
+void
+op_omit::apply(Mat<typename T1::elem_type>& out, const Op<T1,op_omit>& in)
+  {
+  arma_debug_sigprint();
+  
+  typedef typename T1::elem_type eT;
+  
+  const uword omit_mode = in.aux_uword_a;
+  
+  if(arma_config::fast_math_warn)
+    {
+    if(omit_mode == 1)  { arma_warn(1, "omit_nan(): detection of NaN is not reliable in fast math mode"); }
+    if(omit_mode == 2)  { arma_warn(1, "omit_nonfinite(): detection of non-finite values is not reliable in fast math mode"); }
+    }
+  
+  auto is_omitted_1 = [](const eT& x) -> bool { return arma_isnan(x);       };
+  auto is_omitted_2 = [](const eT& x) -> bool { return arma_isnonfinite(x); };
+  
+  if(omit_mode == 1)  { op_omit::apply(out, in.m, is_omitted_1); }
+  if(omit_mode == 2)  { op_omit::apply(out, in.m, is_omitted_2); }
+  }
+
+
+
+template<typename T1, typename functor>
+inline
+void
+op_omit::apply(Mat<typename T1::elem_type>& out, const T1& X, functor is_omitted)
+  {
+  arma_debug_sigprint();
+  
+  typedef typename T1::elem_type eT;
+  
+  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))
+    {
+    const quasi_unwrap<T1> U(X);
+    
+    const eT*   X_mem = U.M.memptr();
+    const uword N     = U.M.n_elem;
+    
+    Mat<eT> Y(N, 1, arma_nozeros_indicator());
+    
+    eT* Y_mem = Y.memptr();
+    
+    uword count = 0;
+    
+    for(uword i=0; i < N; ++i)
+      {
+      const eT val = X_mem[i];
+      
+      if(is_omitted(val) == false)  { Y_mem[count] = val; ++count; }
+      }
+    
+    out.steal_mem_col(Y, count);
+    }
+  else
+    {
+    const Proxy<T1> P(X);
+    
+    const uword N = P.get_n_elem();
+    
+    Mat<eT> Y(N, 1, arma_nozeros_indicator());
+    
+    eT* Y_mem = Y.memptr();
+    
+    uword count = 0;
+    
+    if(Proxy<T1>::use_at == false)
+      {
+      const typename Proxy<T1>::ea_type Pea = P.get_ea();
+      
+      for(uword i=0; i < N; ++i)
+        {
+        const eT val = Pea[i];
+        
+        if(is_omitted(val) == false)  { Y_mem[count] = val; ++count; }
+        }
+      }
+    else
+      {
+      const uword n_rows = P.get_n_rows();
+      const uword n_cols = P.get_n_cols();
+      
+      for(uword c=0; c < n_cols; ++c)
+      for(uword r=0; r < n_rows; ++r)
+        {
+        const eT val = P.at(r,c);
+        
+        if(is_omitted(val) == false)  { Y_mem[count] = val; ++count; }
+        }
+      }
+      
+    out.steal_mem_col(Y, count);
+    }
+  }
+
+
+
+//
+
+
+
+template<typename T1>
+inline
+void
+op_omit_cube::apply(Mat<typename T1::elem_type>& out, const CubeToMatOp<T1, op_omit_cube>& in)
+  {
+  arma_debug_sigprint();
+  
+  typedef typename T1::elem_type eT;
+  
+  const uword omit_mode = in.aux_uword;
+  
+  if(arma_config::fast_math_warn)
+    {
+    if(omit_mode == 1)  { arma_warn(1, "omit_nan(): detection of NaN is not reliable in fast math mode"); }
+    if(omit_mode == 2)  { arma_warn(1, "omit_nonfinite(): detection of non-finite values is not reliable in fast math mode"); }
+    }
+  
+  auto is_omitted_1 = [](const eT& x) -> bool { return arma_isnan(x);       };
+  auto is_omitted_2 = [](const eT& x) -> bool { return arma_isnonfinite(x); };
+  
+  if(omit_mode == 1)  { op_omit_cube::apply(out, in.m, is_omitted_1); }
+  if(omit_mode == 2)  { op_omit_cube::apply(out, in.m, is_omitted_2); }
+  }
+
+
+
+template<typename T1, typename functor>
+inline
+void
+op_omit_cube::apply(Mat<typename T1::elem_type>& out, const T1& X, functor is_omitted)
+  {
+  arma_debug_sigprint();
+  
+  typedef typename T1::elem_type eT;
+  
+  if(is_Cube<T1>::value || is_Cube<typename ProxyCube<T1>::stored_type>::value || (arma_config::openmp && ProxyCube<T1>::use_mp))
+    {
+    const unwrap_cube<T1> U(X);
+    
+    const eT*   X_mem = U.M.memptr();
+    const uword N     = U.M.n_elem;
+    
+    Mat<eT> Y(N, 1, arma_nozeros_indicator());
+    
+    eT* Y_mem = Y.memptr();
+    
+    uword count = 0;
+    
+    for(uword i=0; i < N; ++i)
+      {
+      const eT val = X_mem[i];
+      
+      if(is_omitted(val) == false)  { Y_mem[count] = val; ++count; }
+      }
+    
+    out.steal_mem_col(Y, count);
+    }
+  else
+    {
+    const ProxyCube<T1> P(X);
+    
+    const uword N = P.get_n_elem();
+    
+    Mat<eT> Y(N, 1, arma_nozeros_indicator());
+    
+    eT* Y_mem = Y.memptr();
+    
+    uword count = 0;
+    
+    if(ProxyCube<T1>::use_at == false)
+      {
+      const typename ProxyCube<T1>::ea_type Pea = P.get_ea();
+      
+      for(uword i=0; i < N; ++i)
+        {
+        const eT val = Pea[i];
+        
+        if(is_omitted(val) == false)  { Y_mem[count] = val; ++count; }
+        }
+      }
+    else
+      {
+      const uword n_r = P.get_n_rows();
+      const uword n_c = P.get_n_cols();
+      const uword n_s = P.get_n_slices();
+      
+      for(uword s=0; s < n_s; ++s)
+      for(uword c=0; c < n_c; ++c)
+      for(uword r=0; r < n_r; ++r)
+        {
+        const eT val = P.at(r,c,s);
+        
+        if(is_omitted(val) == false)  { Y_mem[count] = val; ++count; }
+        }
+      }
+      
+    out.steal_mem_col(Y, count);
+    }
+  }
+
+
+
+//! @}
@@ -40,9 +40,18 @@ op_shuffle::apply_direct(Mat<eT>& out, const Mat<eT>& X, const uword dim)
  
  std::vector<packet> packet_vec(N);
  
+  podarray<int> tmp(N);
+  
+  int* tmp_mem = tmp.memptr();
+  
+  const int a = 0;
+  const int b = arma_rng::randi<int>::max_val();
+  
+  arma_rng::randi<int>::fill(tmp_mem, N, a, b);
+  
  for(uword i=0; i<N; ++i)
    {
-    packet_vec[i].val   = int(arma_rng::randi<int>());
+    packet_vec[i].val   = tmp_mem[i];
    packet_vec[i].index = i;
    }
  
@@ -27,10 +27,7 @@ class op_sort_index
  public:
  
  template<typename T1>
-  static inline bool apply_noalias_proxy(Mat<uword>& out, const Proxy<T1>& P, const uword sort_mode);
-  
-  template<typename eT>
-  static inline void apply_noalias_mat(Mat<uword>& out, const Mat<eT>& X, const uword sort_mode);
+  static inline bool apply_helper(Mat<uword>& out, const Proxy<T1>& P, const uword sort_mode);
  
  template<typename T1>
  static inline void apply(Mat<uword>& out, const mtOp<uword,T1,op_sort_index>& in);
@@ -38,20 +35,6 @@ class op_sort_index



-class op_stable_sort_index
-  : public traits_op_col
-  {
-  public:
-  
-  template<typename T1>
-  static inline bool apply_noalias(Mat<uword>& out, const Proxy<T1>& P, const uword sort_mode);
-  
-  template<typename T1>
-  static inline void apply(Mat<uword>& out, const mtOp<uword,T1,op_stable_sort_index>& in);
-  };
-
-
-
 template<typename eT>
 struct arma_sort_index_packet
  {
@@ -98,16 +81,6 @@ struct arma_sort_index_helper_ascend< std::complex<T> >
    {
    return (std::abs(A.val) < std::abs(B.val));
    }
-  
-  // inline
-  // bool
-  // operator() (const arma_sort_index_packet<eT>& A, const arma_sort_index_packet<eT>& B) const
-  //   {
-  //   const T abs_A_val = std::abs(A.val);
-  //   const T abs_B_val = std::abs(B.val);
-  //   
-  //   return ( (abs_A_val != abs_B_val) ? (abs_A_val < abs_B_val) : (std::arg(A.val) < std::arg(B.val)) );
-  //   }
  };


@@ -123,20 +96,14 @@ struct arma_sort_index_helper_descend< std::complex<T> >
    {
    return (std::abs(A.val) > std::abs(B.val));
    }
-  
-  // inline
-  // bool
-  // operator() (const arma_sort_index_packet<eT>& A, const arma_sort_index_packet<eT>& B) const
-  //   {
-  //   const T abs_A_val = std::abs(A.val);
-  //   const T abs_B_val = std::abs(B.val);
-  //   
-  //   return ( (abs_A_val != abs_B_val) ? (abs_A_val > abs_B_val) : (std::arg(A.val) > std::arg(B.val)) );
-  //   }
  };



+//
+
+
+
 template<typename eT>
 struct arma_sort_index_helper_prepare
  {
@@ -21,10 +21,10 @@



-template<typename T1, bool sort_stable>
+template<typename T1>
 inline
 bool
-arma_sort_index_helper(Mat<uword>& out, const Proxy<T1>& P, const uword sort_mode)
+op_sort_index::apply_helper(Mat<uword>& out, const Proxy<T1>& P, const uword sort_mode)
  {
  arma_debug_sigprint();
  
@@ -41,11 +41,13 @@ arma_sort_index_helper(Mat<uword>& out, const Proxy<T1>& P, const uword sort_mod
  
  if(Proxy<T1>::use_at == false)
    {
-    for(uword i=0; i<n_elem; ++i)
+    const typename Proxy<T1>::ea_type Pea = P.get_ea();
+    
+    for(uword i=0; i < n_elem; ++i)
      {
-      const eT val = P[i];
+      const eT val = Pea[i];
      
-      if(arma_isnan(val))  { out.soft_reset(); return false; }
+      if(arma_isnan(val))  { return false; }
      
      packet_vec[i].val   = prepare(val);
      packet_vec[i].index = i;
@@ -63,7 +65,7 @@ arma_sort_index_helper(Mat<uword>& out, const Proxy<T1>& P, const uword sort_mod
      {
      const eT val = P.at(row,col);
      
-      if(arma_isnan(val))  { out.soft_reset(); return false; }
+      if(arma_isnan(val))  { return false; }
      
      packet_vec[i].val   = prepare(val);
      packet_vec[i].index = i;
@@ -79,14 +81,7 @@ arma_sort_index_helper(Mat<uword>& out, const Proxy<T1>& P, const uword sort_mod
    
    arma_sort_index_helper_ascend<T> comparator;
    
-    if(sort_stable == false)
-      {
-      std::sort( packet_vec.begin(), packet_vec.end(), comparator );
-      }
-    else
-      {
-      std::stable_sort( packet_vec.begin(), packet_vec.end(), comparator );
-      }
+    std::stable_sort( packet_vec.begin(), packet_vec.end(), comparator );
    }
  else
    {
@@ -94,14 +89,7 @@ arma_sort_index_helper(Mat<uword>& out, const Proxy<T1>& P, const uword sort_mod
    
    arma_sort_index_helper_descend<T> comparator;
    
-    if(sort_stable == false)
-      {
-      std::sort( packet_vec.begin(), packet_vec.end(), comparator );
-      }
-    else
-      {
-      std::stable_sort( packet_vec.begin(), packet_vec.end(), comparator );
-      }
+    std::stable_sort( packet_vec.begin(), packet_vec.end(), comparator );
    }
  
  uword* out_mem = out.memptr();
@@ -116,40 +104,6 @@ arma_sort_index_helper(Mat<uword>& out, const Proxy<T1>& P, const uword sort_mod



-//
-
-
-
-template<typename T1>
-inline
-bool
-op_sort_index::apply_noalias_proxy(Mat<uword>& out, const Proxy<T1>& P, const uword sort_mode)
-  {
-  arma_debug_sigprint();
-  
-  return arma_sort_index_helper<T1,false>(out, P, sort_mode);
-  }
-
-
-
-template<typename eT>
-inline
-void
-op_sort_index::apply_noalias_mat(Mat<uword>& out, const Mat<eT>& X, const uword sort_mode)
-  {
-  arma_debug_sigprint();
-  
-  if(X.n_elem == 0)  { out.set_size(0,1); return; }
-  
-  const Proxy< Mat<eT> > P(X);
-  
-  const bool all_non_nan = op_sort_index::apply_noalias_proxy(out, P, sort_mode);
-  
-  arma_conform_check( (all_non_nan == false), "sort_index(): detected NaN" );
-  }
-
-
-
 template<typename T1>
 inline
 void
@@ -169,67 +123,20 @@ op_sort_index::apply(Mat<uword>& out, const mtOp<uword,T1,op_sort_index>& in)
    {
    Mat<uword> tmp;
    
-    all_non_nan = op_sort_index::apply_noalias_proxy(tmp, P, sort_mode);
+    all_non_nan = op_sort_index::apply_helper(tmp, P, sort_mode);
    
    out.steal_mem(tmp);
    }
  else
    {
-    all_non_nan = op_sort_index::apply_noalias_proxy(out, P, sort_mode);
+    all_non_nan = op_sort_index::apply_helper(out, P, sort_mode);
    }
  
+  if(all_non_nan == false)  { out.soft_reset(); }
+  
  arma_conform_check( (all_non_nan == false), "sort_index(): detected NaN" );
  }



-//
-
-
-
-template<typename T1>
-inline
-bool
-op_stable_sort_index::apply_noalias(Mat<uword>& out, const Proxy<T1>& P, const uword sort_mode)
-  {
-  arma_debug_sigprint();
-  
-  return arma_sort_index_helper<T1,true>(out, P, sort_mode);
-  }
-
-
-
-template<typename T1>
-inline
-void
-op_stable_sort_index::apply(Mat<uword>& out, const mtOp<uword,T1,op_stable_sort_index>& in)
-  {
-  arma_debug_sigprint();
-  
-  const Proxy<T1> P(in.m);
-  
-  if(P.get_n_elem() == 0)  { out.set_size(0,1); return; }
-  
-  const uword sort_mode = in.aux_uword_a;
-  
-  bool all_non_nan = false;
-  
-  if(P.is_alias(out))
-    {
-    Mat<uword> tmp;
-    
-    all_non_nan = op_stable_sort_index::apply_noalias(tmp, P, sort_mode);
-    
-    out.steal_mem(tmp);
-    }
-  else
-    {
-    all_non_nan = op_stable_sort_index::apply_noalias(out, P, sort_mode);
-    }
-  
-  arma_conform_check( (all_non_nan == false), "stable_sort_index(): detected NaN" );
-  }
-
-
-
 //! @}
@@ -86,7 +86,18 @@ op_sort::apply_noalias(Mat<eT>& out, const Mat<eT>& X, const uword sort_mode, co
        {
        const Col<eT> X_col( const_cast<eT*>(X.colptr(col)), n_rows, false );
        
-        op_sort_index::apply_noalias_mat(indices, X_col, sort_mode);
+        const Proxy< Col<eT> > P(X_col);
+        
+        const bool all_non_nan = op_sort_index::apply_helper(indices, P, sort_mode);
+        
+        if(all_non_nan == false)
+          {
+          out.soft_reset();
+          
+          arma_conform_check(true, "sort(): detected NaN");
+          
+          return;
+          }
        
        const uword* indices_mem = indices.memptr();
        const    eT*   X_col_mem =   X_col.memptr();
@@ -115,7 +126,18 @@ op_sort::apply_noalias(Mat<eT>& out, const Mat<eT>& X, const uword sort_mode, co
        {
        const Col<eT> Y_col( const_cast<eT*>(Y.colptr(col)), n_rows, false );
        
-        op_sort_index::apply_noalias_mat(indices, Y_col, sort_mode);
+        const Proxy< Col<eT> > P(Y_col);
+        
+        const bool all_non_nan = op_sort_index::apply_helper(indices, P, sort_mode);
+        
+        if(all_non_nan == false)
+          {
+          out.soft_reset();
+          
+          arma_conform_check(true, "sort(): detected NaN");
+          
+          return;
+          }
        
        const uword* indices_mem = indices.memptr();
        const    eT*   Y_col_mem =   Y_col.memptr();
@@ -195,9 +217,8 @@ op_sort::apply(Mat<typename T1::elem_type>& out, const Op<T1,op_sort>& in)
  const uword sort_mode = in.aux_uword_a;
  const uword dim       = in.aux_uword_b;
  
-  arma_conform_check( (sort_mode > 1),        "sort(): parameter 'sort_mode' must be 0 or 1" );
-  arma_conform_check( (dim > 1),              "sort(): parameter 'dim' must be 0 or 1"       );
-  arma_conform_check( (X.internal_has_nan()), "sort(): detected NaN"                         );
+  arma_conform_check( (sort_mode > 1), "sort(): parameter 'sort_mode' must be 0 or 1" );
+  arma_conform_check( (dim > 1),       "sort(): parameter 'dim' must be 0 or 1"       );
  
  if(U.is_alias(out))
    {
@@ -229,8 +250,7 @@ op_sort_vec::apply(Mat<typename T1::elem_type>& out, const Op<T1,op_sort_vec>& i
  
  const uword sort_mode = in.aux_uword_a;
  
-  arma_conform_check( (sort_mode > 1),        "sort(): parameter 'sort_mode' must be 0 or 1" );
-  arma_conform_check( (X.internal_has_nan()), "sort(): detected NaN"                         );
+  arma_conform_check( (sort_mode > 1), "sort(): parameter 'sort_mode' must be 0 or 1" );
  
  if(X.n_elem <= 1)  { out = X; return; }
  
@@ -238,7 +258,18 @@ op_sort_vec::apply(Mat<typename T1::elem_type>& out, const Op<T1,op_sort_vec>& i
    {
    uvec indices;
    
-    op_sort_index::apply_noalias_mat(indices, X, sort_mode);
+    const Proxy< Mat<eT> > P(X);
+    
+    const bool all_non_nan = op_sort_index::apply_helper(indices, P, sort_mode);
+    
+    if(all_non_nan == false)
+      {
+      out.soft_reset();
+      
+      arma_conform_check(true, "sort(): detected NaN");
+      
+      return;
+      }
    
    const uword N = indices.n_elem;
    
@@ -31,6 +31,37 @@ op_sp_nonzeros::apply(Mat<typename T1::elem_type>& out, const SpToDOp<T1, op_sp_
  
  typedef typename T1::elem_type eT;
  
+  if(is_SpMat<T1>::value || is_SpMat<typename SpProxy<T1>::stored_type>::value)
+    {
+    const unwrap_spmat<T1> U(X.m);
+    
+    out.set_size(U.M.n_nonzero,1);
+    
+    arrayops::copy(out.memptr(), U.M.values, U.M.n_nonzero);
+    
+    return;
+    }
+  
+  if(is_SpSubview<T1>::value)
+    {
+    const SpSubview<eT>& sv = reinterpret_cast< const SpSubview<eT>& >(X.m);
+    
+    if(sv.n_rows == sv.m.n_rows)
+      {
+      arma_debug_print("op_sp_nonzeros::apply(): SpSubview optimisation");
+      
+      const SpMat<eT>& m   = sv.m;
+      const uword      col = sv.aux_col1;
+      const uword      N   = sv.n_nonzero;
+      
+      out.set_size(N, 1);
+      
+      arrayops::copy(out.memptr(), &(m.values[ m.col_ptrs[col] ]), N);
+      
+      return;
+      }
+    }
+  
  const SpProxy<T1> P(X.m);
  
  const uword N = P.get_n_nonzero();
@@ -39,30 +70,6 @@ op_sp_nonzeros::apply(Mat<typename T1::elem_type>& out, const SpToDOp<T1, op_sp_
  
  if(N == 0)  { return; }
  
-  if(is_SpMat<typename SpProxy<T1>::stored_type>::value)
-    {
-    const unwrap_spmat<typename SpProxy<T1>::stored_type> U(P.Q);
-    
-    arrayops::copy(out.memptr(), U.M.values, N);
-    
-    return;
-    }
-  
-  if(is_SpSubview<typename SpProxy<T1>::stored_type>::value)
-    {
-    const SpSubview<eT>& sv = reinterpret_cast< const SpSubview<eT>& >(P.Q);
-    
-    if(sv.n_rows == sv.m.n_rows)
-      {
-      const SpMat<eT>& m   = sv.m;
-      const uword      col = sv.aux_col1;
-      
-      arrayops::copy(out.memptr(), &(m.values[ m.col_ptrs[col] ]), N);
-      
-      return;
-      }
-    }
-  
  eT* out_mem = out.memptr();
  
  typename SpProxy<T1>::const_iterator_type it = P.begin();
@@ -25,11 +25,11 @@ class op_sp_sum
  {
  public:
  
-  template<typename T1>
-  inline static void apply(Mat<typename T1::elem_type>& out, const mtSpReduceOp<typename T1::elem_type, T1, op_sp_sum>& in);
+  template<typename eT, typename T1>
+  inline static void apply(Mat<eT>& out, const mtSpReduceOp<eT, T1, op_sp_sum>& in);
  
-  template<typename T1>
-  inline static void apply(Mat<typename T1::elem_type>& out, const mtSpReduceOp<typename T1::elem_type, SpOp<T1, spop_square>, op_sp_sum>& in);
+  template<typename eT, typename T1>
+  inline static void apply(Mat<eT>& out, const mtSpReduceOp<eT, SpOp<T1, spop_square>, op_sp_sum>& in);
  };


@@ -21,14 +21,14 @@



-template<typename T1>
+template<typename eT, typename T1>
 inline
 void
-op_sp_sum::apply(Mat<typename T1::elem_type>& out, const mtSpReduceOp<typename T1::elem_type, T1, op_sp_sum>& in)
+op_sp_sum::apply(Mat<eT>& out, const mtSpReduceOp<eT, T1, op_sp_sum>& in)
  {
  arma_debug_sigprint();
  
-  typedef typename T1::elem_type eT;
+  // NOTE: assuming that eT is the same as T1::elem_type
  
  const uword dim = in.aux_uword_a;
  
@@ -83,14 +83,14 @@ op_sp_sum::apply(Mat<typename T1::elem_type>& out, const mtSpReduceOp<typename T



-template<typename T1>
+template<typename eT, typename T1>
 inline
 void
-op_sp_sum::apply(Mat<typename T1::elem_type>& out, const mtSpReduceOp<typename T1::elem_type, SpOp<T1, spop_square>, op_sp_sum>& in)
+op_sp_sum::apply(Mat<eT>& out, const mtSpReduceOp<eT, SpOp<T1, spop_square>, op_sp_sum>& in)
  {
  arma_debug_sigprint();
  
-  typedef typename T1::elem_type eT;
+  // NOTE: assuming that eT is the same as T1::elem_type
  
  const uword dim = in.aux_uword_a;
  
@@ -32,16 +32,8 @@ class op_stddev
  template<typename in_eT>
  inline static void apply_noalias(Mat<typename get_pod_type<in_eT>::result>& out, const Mat<in_eT>& X, const uword norm_type, const uword dim);
  
-  //
-  
-  template<typename eT>
-  inline static typename get_pod_type<eT>::result stddev_vec(const subview_col<eT>& X, const uword norm_type = 0);
-  
-  template<typename eT>
-  inline static typename get_pod_type<eT>::result stddev_vec(const subview_row<eT>& X, const uword norm_type = 0);
-  
  template<typename T1>
-  inline static typename T1::pod_type stddev_vec(const Base<typename T1::elem_type, T1>& X, const uword norm_type = 0);
+  inline static typename T1::pod_type stddev_vec(const T1& X, const uword norm_type = 0);
  };


@@ -76,7 +76,7 @@ op_stddev::apply_noalias(Mat<typename get_pod_type<in_eT>::result>& out, const M
      {
      out_eT* out_mem = out.memptr();
      
-      for(uword col=0; col<X_n_cols; ++col)
+      for(uword col=0; col < X_n_cols; ++col)
        {
        out_mem[col] = std::sqrt( op_var::direct_var( X.colptr(col), X_n_rows, norm_type ) );
        }
@@ -91,16 +91,15 @@ op_stddev::apply_noalias(Mat<typename get_pod_type<in_eT>::result>& out, const M
    
    if(X_n_cols > 0)
      {
-      podarray<in_eT> dat(X_n_cols);
-      
-      in_eT*  dat_mem = dat.memptr();
      out_eT* out_mem = out.memptr();
      
-      for(uword row=0; row<X_n_rows; ++row)
+      podarray<in_eT> tmp;
+      
+      for(uword row=0; row < X_n_rows; ++row)
        {
-        dat.copy_row(X, row);
+        tmp.copy_row(X, row);
        
-        out_mem[row] = std::sqrt( op_var::direct_var( dat_mem, X_n_cols, norm_type) );
+        out_mem[row] = std::sqrt( op_var::direct_var( tmp.memptr(), tmp.n_elem, norm_type) );
        }
      }
    }
@@ -111,7 +110,7 @@ op_stddev::apply_noalias(Mat<typename get_pod_type<in_eT>::result>& out, const M
 template<typename T1>
 inline
 typename T1::pod_type
-op_stddev::stddev_vec(const Base<typename T1::elem_type, T1>& X, const uword norm_type)
+op_stddev::stddev_vec(const T1& X, const uword norm_type)
  {
  arma_debug_sigprint();
  
@@ -119,7 +118,7 @@ op_stddev::stddev_vec(const Base<typename T1::elem_type, T1>& X, const uword nor
  
  arma_conform_check( (norm_type > 1), "stddev(): parameter 'norm_type' must be 0 or 1" );
  
-  const quasi_unwrap<T1> U(X.get_ref());
+  const quasi_unwrap<T1> U(X);
  
  if(U.M.n_elem == 0)
    {
@@ -133,67 +132,4 @@ op_stddev::stddev_vec(const Base<typename T1::elem_type, T1>& X, const uword nor



-template<typename eT>
-inline
-typename get_pod_type<eT>::result
-op_stddev::stddev_vec(const subview_col<eT>& X, const uword norm_type)
-  {
-  arma_debug_sigprint();
-  
-  typedef typename get_pod_type<eT>::result T;
-  
-  arma_conform_check( (norm_type > 1), "stddev(): parameter 'norm_type' must be 0 or 1" );
-  
-  if(X.n_elem == 0)
-    {
-    arma_conform_check(true, "stddev(): object has no elements");
-    
-    return Datum<T>::nan;
-    }
-  
-  return std::sqrt( op_var::direct_var(X.colptr(0), X.n_rows, norm_type) );
-  }
-
-
-
-
-template<typename eT>
-inline
-typename get_pod_type<eT>::result
-op_stddev::stddev_vec(const subview_row<eT>& X, const uword norm_type)
-  {
-  arma_debug_sigprint();
-  
-  typedef typename get_pod_type<eT>::result T;
-  
-  arma_conform_check( (norm_type > 1), "stddev(): parameter 'norm_type' must be 0 or 1" );
-  
-  if(X.n_elem == 0)
-    {
-    arma_conform_check(true, "stddev(): object has no elements");
-    
-    return Datum<T>::nan;
-    }
-  
-  const Mat<eT>& A = X.m;
-  
-  const uword start_row = X.aux_row1;
-  const uword start_col = X.aux_col1;
-  
-  const uword end_col_p1 = start_col + X.n_cols;
-  
-  podarray<eT> tmp(X.n_elem);
-  eT* tmp_mem = tmp.memptr();
-  
-  for(uword i=0, col=start_col; col < end_col_p1; ++col, ++i)
-    {
-    tmp_mem[i] = A.at(start_row, col);
-    }
-  
-  return std::sqrt( op_var::direct_var(tmp.memptr(), tmp.n_elem, norm_type) );
-  }
-
-
-
 //! @}
-
@@ -32,10 +32,13 @@ class op_sum
  
  template<typename T1>
  inline static void apply(Mat<typename T1::elem_type>& out, const Op< eOp<T1,eop_square>, op_sum >& in);
-
+  
  template<typename T1>
  inline static void apply(Mat<typename T1::elem_type>& out, const Op< eOp<T1,eop_pow   >, op_sum >& in);
-
+  
+  template<typename T1>
+  inline static void apply_generic(Mat<typename T1::elem_type>& out, const Op<T1,op_sum>& in);
+  
  template<typename eT>
  inline static void apply_mat_noalias(Mat<eT>& out, const Mat<eT>& X, const uword dim);
  
@@ -28,6 +28,92 @@ op_sum::apply(Mat<typename T1::elem_type>& out, const Op<T1,op_sum>& in)
  {
  arma_debug_sigprint();
  
+  op_sum::apply_generic(out, in);
+  }
+
+
+
+template<typename T1>
+inline
+void
+op_sum::apply(Mat<typename T1::elem_type>& out, const Op< eOp<T1,eop_square>, op_sum >& in)
+  {
+  arma_debug_sigprint();
+  
+  typedef typename T1::elem_type eT;
+  
+  typedef eOp<T1,eop_square> inner_expr_type;
+  
+  typedef typename inner_expr_type::proxy_type::stored_type inner_expr_P_stored_type;
+  
+  if(is_Mat<inner_expr_P_stored_type>::value)
+    {
+    const uword dim = in.aux_uword_a;
+    
+    arma_conform_check( (dim > 1), "sum(): parameter 'dim' must be 0 or 1" );
+    
+    const quasi_unwrap<inner_expr_P_stored_type> U(in.m.P.Q);
+    
+    if(U.is_alias(out))
+      {
+      Mat<eT> tmp;
+      
+      op_sum::apply_mat_square_noalias(tmp, U.M, dim);
+      
+      out.steal_mem(tmp);
+      }
+    else
+      {
+      op_sum::apply_mat_square_noalias(out, U.M, dim);
+      }
+    
+    return;
+    }
+  
+  op_sum::apply_generic(out, in);
+  }
+
+
+
+template<typename T1>
+inline
+void
+op_sum::apply(Mat<typename T1::elem_type>& out, const Op< eOp<T1,eop_pow>, op_sum >& in)
+  {
+  arma_debug_sigprint();
+  
+  typedef typename T1::elem_type eT;
+  
+  if(arma_config::optimise_powexpr && (in.m.aux == eT(2)))
+    {
+    typedef Op< eOp<T1,eop_square>, op_sum > modified_whole_expr_type;
+    
+    op_sum::apply(out, reinterpret_cast<const modified_whole_expr_type& >(in) );
+    
+    return;
+    }
+  
+  if(arma_config::optimise_powexpr && (in.m.aux == eT(0.5)) && is_non_integral<eT>::value)
+    {
+    typedef Op< eOp<T1,eop_sqrt>, op_sum > modified_whole_expr_type;
+    
+    op_sum::apply(out, reinterpret_cast<const modified_whole_expr_type& >(in) );
+    
+    return;
+    }
+  
+  op_sum::apply_generic(out, in);
+  }
+
+
+
+template<typename T1>
+inline
+void
+op_sum::apply_generic(Mat<typename T1::elem_type>& out, const Op<T1,op_sum>& in)
+  {
+  arma_debug_sigprint();
+  
  typedef typename T1::elem_type eT;
  
  const uword dim = in.aux_uword_a;
@@ -72,152 +158,6 @@ op_sum::apply(Mat<typename T1::elem_type>& out, const Op<T1,op_sum>& in)



-template<typename T1>
-inline
-void
-op_sum::apply(Mat<typename T1::elem_type>& out, const Op< eOp<T1,eop_square>, op_sum >& in)
-  {
-  arma_debug_sigprint();
-  
-  typedef typename T1::elem_type eT;
-  
-  typedef eOp<T1,eop_square> inner_expr_type;
-  
-  typedef typename inner_expr_type::proxy_type::stored_type inner_expr_P_stored_type;
-  
-  const uword dim = in.aux_uword_a;
-  
-  arma_conform_check( (dim > 1), "sum(): parameter 'dim' must be 0 or 1" );
-  
-  if(is_Mat<inner_expr_P_stored_type>::value)
-    {
-    const quasi_unwrap<inner_expr_P_stored_type> U(in.m.P.Q);
-    
-    if(U.is_alias(out))
-      {
-      Mat<eT> tmp;
-      
-      op_sum::apply_mat_square_noalias(tmp, U.M, dim);
-      
-      out.steal_mem(tmp);
-      }
-    else
-      {
-      op_sum::apply_mat_square_noalias(out, U.M, dim);
-      }
-    }
-  else
-  if(arma_config::openmp && Proxy<inner_expr_type>::use_mp)
-    {
-    const quasi_unwrap<inner_expr_type> U(in.m);  // force evaluation of compound inner expression
-    
-    if(U.is_alias(out))
-      {
-      Mat<eT> tmp;
-      
-      op_sum::apply_mat_noalias(tmp, U.M, dim);
-      
-      out.steal_mem(tmp);
-      }
-    else
-      {
-      op_sum::apply_mat_noalias(out, U.M, dim);
-      }
-    }
-  else
-    {
-    const Proxy<inner_expr_type> P(in.m);
-    
-    if(P.is_alias(out))
-      {
-      Mat<eT> tmp;
-      
-      op_sum::apply_proxy_noalias(tmp, P, dim);
-      
-      out.steal_mem(tmp);
-      }
-    else
-      {
-      op_sum::apply_proxy_noalias(out, P, dim);
-      }
-    }
-  }
-
-
-
-template<typename T1>
-inline
-void
-op_sum::apply(Mat<typename T1::elem_type>& out, const Op< eOp<T1,eop_pow>, op_sum >& in)
-  {
-  arma_debug_sigprint();
-  
-  typedef typename T1::elem_type eT;
-  
-  if(in.m.aux == eT(2))
-    {
-    typedef Op< eOp<T1,eop_square>, op_sum > modified_whole_expr_type;
-    
-    op_sum::apply(out, reinterpret_cast<const modified_whole_expr_type& >(in) );
-    
-    return;
-    }
-  
-  if((in.m.aux == eT(0.5)) && is_non_integral<eT>::value)
-    {
-    typedef Op< eOp<T1,eop_sqrt>, op_sum > modified_whole_expr_type;
-    
-    op_sum::apply(out, reinterpret_cast<const modified_whole_expr_type& >(in) );
-    
-    return;
-    }
-  
-  typedef eOp<T1,eop_pow> inner_expr_type;
-  
-  typedef typename inner_expr_type::proxy_type::stored_type inner_expr_P_stored_type;
-  
-  const uword dim = in.aux_uword_a;
-  
-  arma_conform_check( (dim > 1), "sum(): parameter 'dim' must be 0 or 1" );
-  
-  if( (is_Mat<inner_expr_P_stored_type>::value) || (arma_config::openmp && Proxy<inner_expr_type>::use_mp) )
-    {
-    const quasi_unwrap<inner_expr_type> U(in.m);  // force evaluation of eop_pow
-    
-    if(U.is_alias(out))
-      {
-      Mat<eT> tmp;
-      
-      op_sum::apply_mat_noalias(tmp, U.M, dim);
-      
-      out.steal_mem(tmp);
-      }
-    else
-      {
-      op_sum::apply_mat_noalias(out, U.M, dim);
-      }
-    }
-  else
-    {
-    const Proxy<inner_expr_type> P(in.m);
-    
-    if(P.is_alias(out))
-      {
-      Mat<eT> tmp;
-      
-      op_sum::apply_proxy_noalias(tmp, P, dim);
-      
-      out.steal_mem(tmp);
-      }
-    else
-      {
-      op_sum::apply_proxy_noalias(out, P, dim);
-      }
-    }
-  }
-
-
-
 template<typename eT>
 inline
 void
@@ -34,23 +34,17 @@ class op_var
  
  //
  
-  template<typename eT>
-  inline static typename get_pod_type<eT>::result var_vec(const subview_col<eT>& X, const uword norm_type = 0);
-  
-  template<typename eT>
-  inline static typename get_pod_type<eT>::result var_vec(const subview_row<eT>& X, const uword norm_type = 0);
-  
  template<typename T1>
-  inline static typename T1::pod_type var_vec(const Base<typename T1::elem_type, T1>& X, const uword norm_type = 0);
+  inline static typename T1::pod_type var_vec(const T1& X, const uword norm_type = 0);
  
  
  //
  
  template<typename eT>
-  inline static eT direct_var(const eT* const X, const uword N, const uword norm_type = 0);
+  inline static eT direct_var(const eT* X, const uword N, const uword norm_type = 0);
  
  template<typename eT>
-  inline static eT direct_var_robust(const eT* const X, const uword N, const uword norm_type = 0);
+  inline static eT direct_var_robust(const eT* X, const uword N, const uword norm_type = 0);
  
  
  //
@@ -76,7 +76,7 @@ op_var::apply_noalias(Mat<typename get_pod_type<in_eT>::result>& out, const Mat<
      {
      out_eT* out_mem = out.memptr();
      
-      for(uword col=0; col<X_n_cols; ++col)
+      for(uword col=0; col < X_n_cols; ++col)
        {
        out_mem[col] = op_var::direct_var( X.colptr(col), X_n_rows, norm_type );
        }
@@ -91,16 +91,15 @@ op_var::apply_noalias(Mat<typename get_pod_type<in_eT>::result>& out, const Mat<
    
    if(X_n_cols > 0)
      {
-      podarray<in_eT> dat(X_n_cols);
-      
-      in_eT*  dat_mem = dat.memptr();
      out_eT* out_mem = out.memptr();
      
-      for(uword row=0; row<X_n_rows; ++row)
+      podarray<in_eT> tmp;
+      
+      for(uword row=0; row < X_n_rows; ++row)
        {
-        dat.copy_row(X, row);
+        tmp.copy_row(X, row);
        
-        out_mem[row] = op_var::direct_var( dat_mem, X_n_cols, norm_type );
+        out_mem[row] = op_var::direct_var( tmp.memptr(), tmp.n_elem, norm_type );
        }
      }
    }
@@ -111,7 +110,7 @@ op_var::apply_noalias(Mat<typename get_pod_type<in_eT>::result>& out, const Mat<
 template<typename T1>
 inline
 typename T1::pod_type
-op_var::var_vec(const Base<typename T1::elem_type, T1>& X, const uword norm_type)
+op_var::var_vec(const T1& X, const uword norm_type)
  {
  arma_debug_sigprint();
  
@@ -119,7 +118,7 @@ op_var::var_vec(const Base<typename T1::elem_type, T1>& X, const uword norm_type
  
  arma_conform_check( (norm_type > 1), "var(): parameter 'norm_type' must be 0 or 1" );
  
-  const quasi_unwrap<T1> U(X.get_ref());
+  const quasi_unwrap<T1> U(X);
  
  if(U.M.n_elem == 0)
    {
@@ -133,73 +132,11 @@ op_var::var_vec(const Base<typename T1::elem_type, T1>& X, const uword norm_type



-template<typename eT>
-inline
-typename get_pod_type<eT>::result
-op_var::var_vec(const subview_col<eT>& X, const uword norm_type)
-  {
-  arma_debug_sigprint();
-  
-  typedef typename get_pod_type<eT>::result T;
-  
-  arma_conform_check( (norm_type > 1), "var(): parameter 'norm_type' must be 0 or 1" );
-  
-  if(X.n_elem == 0)
-    {
-    arma_conform_check(true, "var(): object has no elements");
-    
-    return Datum<T>::nan;
-    }
-  
-  return op_var::direct_var(X.colptr(0), X.n_rows, norm_type);
-  }
-
-
-
-
-template<typename eT>
-inline
-typename get_pod_type<eT>::result
-op_var::var_vec(const subview_row<eT>& X, const uword norm_type)
-  {
-  arma_debug_sigprint();
-  
-  typedef typename get_pod_type<eT>::result T;
-  
-  arma_conform_check( (norm_type > 1), "var(): parameter 'norm_type' must be 0 or 1" );
-  
-  if(X.n_elem == 0)
-    {
-    arma_conform_check(true, "var(): object has no elements");
-    
-    return Datum<T>::nan;
-    }
-  
-  const Mat<eT>& A = X.m;
-  
-  const uword start_row = X.aux_row1;
-  const uword start_col = X.aux_col1;
-  
-  const uword end_col_p1 = start_col + X.n_cols;
-  
-  podarray<eT> tmp(X.n_elem);
-  eT* tmp_mem = tmp.memptr();
-  
-  for(uword i=0, col=start_col; col < end_col_p1; ++col, ++i)
-    {
-    tmp_mem[i] = A.at(start_row, col);
-    }
-  
-  return op_var::direct_var(tmp.memptr(), tmp.n_elem, norm_type);
-  }
-
-
-
 //! find the variance of an array
 template<typename eT>
 inline
 eT
-op_var::direct_var(const eT* const X, const uword n_elem, const uword norm_type)
+op_var::direct_var(const eT* X, const uword n_elem, const uword norm_type)
  {
  arma_debug_sigprint();
  
@@ -207,6 +144,8 @@ op_var::direct_var(const eT* const X, const uword n_elem, const uword norm_type)
    {
    const eT acc1 = op_mean::direct_mean(X, n_elem);
    
+    if(arma_isnonfinite(acc1))  { return Datum<eT>::nan; }
+    
    eT acc2 = eT(0);
    eT acc3 = eT(0);
    
@@ -251,7 +190,7 @@ op_var::direct_var(const eT* const X, const uword n_elem, const uword norm_type)
 template<typename eT>
 inline
 eT
-op_var::direct_var_robust(const eT* const X, const uword n_elem, const uword norm_type)
+op_var::direct_var_robust(const eT* X, const uword n_elem, const uword norm_type)
  {
  arma_debug_sigprint();
  
@@ -294,6 +233,8 @@ op_var::direct_var(const std::complex<T>* const X, const uword n_elem, const uwo
    {
    const eT acc1 = op_mean::direct_mean(X, n_elem);
    
+    if(arma_isnonfinite(acc1))  { return Datum<T>::nan; }
+    
    T  acc2 =  T(0);
    eT acc3 = eT(0);
    
@@ -354,4 +295,3 @@ op_var::direct_var_robust(const std::complex<T>* const X, const uword n_elem, co


 //! @}
-
@@ -37,8 +37,8 @@ podarray<eT>::~podarray()
 template<typename eT>
 inline
 podarray<eT>::podarray()
-  : n_elem(0)
-  , mem   (0)
+  : n_elem(0      )
+  , mem   (nullptr)
  {
  arma_debug_sigprint_this(this);
  }
@@ -128,7 +128,7 @@ arma_inline
 eT&
 podarray<eT>::operator[] (const uword i)
  {
-  return access::rw(mem[i]);
+  return mem[i];
  }


@@ -152,7 +152,7 @@ podarray<eT>::operator() (const uword i)
  {
  arma_conform_check_bounds( (i >= n_elem), "podarray::operator(): index out of bounds" );
  
-  return access::rw(mem[i]);
+  return mem[i];
  }


@@ -258,11 +258,11 @@ podarray<eT>::copy_row(const Mat<eT>& A, const uword row)
  {
  arma_debug_sigprint();
  
-  // note: this function assumes that the podarray has been set to the correct size beforehand
-  
  const uword n_rows = A.n_rows;
  const uword n_cols = A.n_cols;
  
+  init_warm(n_cols);
+  
  const eT*   A_mem = &(A.at(row,0));
        eT* out_mem = memptr();
  
@@ -162,9 +162,9 @@ running_stat<eT>::operator() (const typename running_stat<eT>::T sample)
  {
  arma_debug_sigprint();
  
-  if( arma_isfinite(sample) == false )
+  if(arma_isnonfinite(sample))
    {
-    arma_warn(3, "running_stat: sample ignored as it is non-finite" );
+    arma_warn(3, "running_stat: non-finite sample ignored" );
    return;
    }
  
@@ -181,9 +181,9 @@ running_stat<eT>::operator() (const std::complex< typename running_stat<eT>::T >
  {
  arma_debug_sigprint();
  
-  if( arma_isfinite(sample) == false )
+  if(arma_isnonfinite(sample))
    {
-    arma_warn(3, "running_stat: sample ignored as it is non-finite" );
+    arma_warn(3, "running_stat: non-finite sample ignored" );
    return;
    }
  
@@ -100,7 +100,7 @@ running_stat_vec<obj_type>::operator() (const Base<typename running_stat_vec<obj
  
  if( sample.internal_has_nonfinite() )
    {
-    arma_warn(3, "running_stat_vec: sample ignored as it has non-finite elements");
+    arma_warn(3, "running_stat_vec: non-finite sample ignored");
    return;
    }
  
@@ -128,7 +128,7 @@ running_stat_vec<obj_type>::operator() (const Base< std::complex<typename runnin
  
  if( sample.internal_has_nonfinite() )
    {
-    arma_warn(3, "running_stat_vec: sample ignored as it has non-finite elements");
+    arma_warn(3, "running_stat_vec: non-finite sample ignored");
    return;
    }
  
@@ -1237,21 +1237,21 @@ sp_auxlib::spsolve_simple(Mat<typename T1::elem_type>& X, const SpBase<typename
    
    superlu_stat_wrangler stat;
    
-    int info = 0; // Return code.
+    superlu::int_t info = 0; // Return code.
    
    arma_debug_print("superlu::gssv()");
    superlu::gssv<eT>(&options, a.get_ptr(), perm_c.get_ptr(), perm_r.get_ptr(), l.get_ptr(), u.get_ptr(), x.get_ptr(), stat.get_ptr(), &info);
    
    
    // Process the return code.
-    if( (info > 0) && (info <= int(A.n_cols)) )
+    if( (info > 0) && (info <= superlu::int_t(A.n_cols)) )
      {
      // std::ostringstream tmp;
      // tmp << "spsolve(): could not solve system; LU factorisation completed, but detected zero in U(" << (info-1) << ',' << (info-1) << ')';
      // arma_warn(1, tmp.str());
      }
    else
-    if(info > int(A.n_cols))
+    if(info > superlu::int_t(A.n_cols))
      {
      arma_warn(1, "spsolve(): memory allocation failure");
      }
@@ -1381,10 +1381,11 @@ sp_auxlib::spsolve_refine(Mat<typename T1::elem_type>& X, typename T1::pod_type&
    char equed[8] = {};     // extra characters for paranoia
    T    rpg      = T(0);
    T    rcond    = T(0);
-    int  info     = int(0); // Return code.
    
-    char  work[8] = {};
-    int  lwork    = int(0);  // 0 means superlu will allocate memory
+    char work[8] = {};
+    superlu::int_t lwork = 0;  // 0 means superlu will allocate memory
+    
+    superlu::int_t info = 0; // Return code.
    
    arma_debug_print("superlu::gssvx()");
    superlu::gssvx<eT>(&options, a.get_ptr(), perm_c.get_ptr(), perm_r.get_ptr(), etree.get_ptr(), equed, R.get_ptr(), C.get_ptr(), l.get_ptr(), u.get_ptr(), &work[0], lwork, b.get_ptr(), x.get_ptr(), &rpg, &rcond, ferr.get_ptr(), berr.get_ptr(), &glu, &mu, stat.get_ptr(), &info);
@@ -1396,20 +1397,20 @@ sp_auxlib::spsolve_refine(Mat<typename T1::elem_type>& X, typename T1::pod_type&
      {
      status = true;
      }
-    if( (info > 0) && (info <= int(A.n_cols)) )
+    if( (info > 0) && (info <= superlu::int_t(A.n_cols)) )
      {
      // std::ostringstream tmp;
      // tmp << "spsolve(): could not solve system; LU factorisation completed, but detected zero in U(" << (info-1) << ',' << (info-1) << ')';
      // arma_warn(1, tmp.str());
      }
    else
-    if( (info == int(A.n_cols+1)) && (user_opts.allow_ugly) )
+    if( (info == superlu::int_t(A.n_cols+1)) && (user_opts.allow_ugly) )
      {
      arma_warn(2, "spsolve(): system is singular to working precision (rcond: ", rcond, ")");
      status = true;
      }
    else
-    if(info > int(A.n_cols+1))
+    if(info > superlu::int_t(A.n_cols+1))
      {
      arma_warn(1, "spsolve(): memory allocation failure");
      }
@@ -2027,6 +2028,7 @@ sp_auxlib::run_aupd_plain
    arma_ignore(n_eigvals);
    arma_ignore(which);
    arma_ignore(X);
+    arma_ignore(Xst);
    arma_ignore(sym);
    arma_ignore(n);
    arma_ignore(tol);
@@ -2123,7 +2125,7 @@ sp_auxlib::run_aupd_shiftinvert
    superlu_opts superlu_opts_default;
    superlu::superlu_options_t options;
    sp_auxlib::set_superlu_opts(options, superlu_opts_default);
-    int lwork = 0;
+    
    superlu::trans_t trans = superlu::NOTRANS;
    
    superlu::GlobalLU_t Glu; /* Not needed on return. */
@@ -2176,7 +2178,9 @@ sp_auxlib::run_aupd_shiftinvert
    
    int panel_size = superlu::sp_ispec_environ(1);
    int relax      = superlu::sp_ispec_environ(2);
-    int slu_info   = 0; // Return code.
+    
+    superlu::int_t lwork    = 0;
+    superlu::int_t slu_info = 0; // Return code.
    
    arma_debug_print("superlu::gstrf()");
    superlu::get_permutation_c(options.ColPerm, x.get_ptr(), perm_c.get_ptr());
@@ -2704,8 +2708,9 @@ superlu_worker<eT>::factorise(typename get_pod_type<eT>::result& out_rcond, cons
  
  int panel_size = superlu::sp_ispec_environ(1);
  int relax      = superlu::sp_ispec_environ(2);
-  int lwork      = 0;
-  int info       = 0;
+  
+  superlu::int_t lwork = 0;
+  superlu::int_t info  = 0;
  
  arma_debug_print("superlu::superlu::get_permutation_c()");
  superlu::get_permutation_c(options.ColPerm, AA.get_ptr(), perm_c.get_ptr());
@@ -134,8 +134,6 @@ spdiagview<eT>::operator*=(const eT val)
  {
  arma_debug_sigprint();
  
-  if(val == eT(0))  { (*this).zeros(); return; }
-  
  SpMat<eT>& t_m = const_cast< SpMat<eT>& >(m);
  
  const uword t_n_elem     = n_elem;
@@ -60,13 +60,13 @@ spglue_schur::apply_noalias(SpMat<eT>& out, const SpProxy<T1>& pa, const SpProxy
  
  arma_conform_assert_same_size(pa.get_n_rows(), pa.get_n_cols(), pb.get_n_rows(), pb.get_n_cols(), "element-wise multiplication");
  
-  if( (pa.get_n_nonzero() == 0) || (pb.get_n_nonzero() == 0) )
+  if( (pa.get_n_nonzero() == 0) && (pb.get_n_nonzero() == 0) )
    {
    out.zeros(pa.get_n_rows(), pa.get_n_cols());
    return;
    }
  
-  const uword max_n_nonzero = (std::min)(pa.get_n_nonzero(), pb.get_n_nonzero());
+  const uword max_n_nonzero = pa.get_n_nonzero() + pb.get_n_nonzero();
  
  // Resize memory to upper bound
  out.reserve(pa.get_n_rows(), pa.get_n_cols(), max_n_nonzero);
@@ -82,24 +82,19 @@ spglue_schur::apply_noalias(SpMat<eT>& out, const SpProxy<T1>& pa, const SpProxy
  
  while( (x_it != x_end) || (y_it != y_end) )
    {
+    eT out_val;
+    
    const uword x_it_row = x_it.row();
    const uword x_it_col = x_it.col();
    
    const uword y_it_row = y_it.row();
    const uword y_it_col = y_it.col();
    
+    bool use_y_loc = false;
+    
    if(x_it == y_it)
      {
-      const eT out_val = (*x_it) * (*y_it);
-      
-      if(out_val != eT(0))
-        {
-        access::rw(out.values[count]) = out_val;
-        
-        access::rw(out.row_indices[count]) = x_it_row;
-        access::rw(out.col_ptrs[x_it_col + 1])++;
-        ++count;
-        }
+      out_val = (*x_it) * (*y_it);
      
      ++x_it;
      ++y_it;
@@ -108,14 +103,32 @@ spglue_schur::apply_noalias(SpMat<eT>& out, const SpProxy<T1>& pa, const SpProxy
      {
      if((x_it_col < y_it_col) || ((x_it_col == y_it_col) && (x_it_row < y_it_row))) // if y is closer to the end
        {
+        out_val = (*x_it) * eT(0);  // in case (*x_it) is inf or nan
+        
        ++x_it;
        }
      else
        {
+        out_val = eT(0) * (*y_it);  // in case (*y_it) is inf or nan
+        
        ++y_it;
+        
+        use_y_loc = true;
        }
      }
    
+    if(out_val != eT(0))
+      {
+      access::rw(out.values[count]) = out_val;
+      
+      const uword out_row = (use_y_loc == false) ? x_it_row : y_it_row;
+      const uword out_col = (use_y_loc == false) ? x_it_col : y_it_col;
+      
+      access::rw(out.row_indices[count]) = out_row;
+      access::rw(out.col_ptrs[out_col + 1])++;
+      ++count;
+      }
+    
    arma_check( (count > max_n_nonzero), "internal error: spglue_schur::apply_noalias(): count > max_n_nonzero" );
    }
  
@@ -262,4 +262,15 @@ class spop_fliplr



+class spop_replace
+  : public traits_op_passthru
+  {
+  public:
+  
+  template<typename eT, typename T1>
+  inline static void apply(SpMat<eT>& out, const mtSpOp<eT, T1, spop_replace>& in);
+  };
+
+
+
 //! @}
@@ -45,16 +45,7 @@ spop_scalar_times::apply(SpMat<typename T1::elem_type>& out, const SpOp<T1,spop_
  
  typedef typename T1::elem_type eT;
  
-  if(in.aux != eT(0))
-    {
-    out.init_xform(in.m, priv::functor_scalar_times<eT>(in.aux));
-    }
-  else
-    {
-    const SpProxy<T1> P(in.m);
-    
-    out.zeros( P.get_n_rows(), P.get_n_cols() );
-    }
+  out.init_xform(in.m, priv::functor_scalar_times<eT>(in.aux));
  }


@@ -83,19 +74,9 @@ spop_cx_scalar_times::apply(SpMat< std::complex<typename T1::pod_type> >& out, c
  {
  arma_debug_sigprint();
  
-  typedef typename T1::pod_type         T;
-  typedef typename std::complex<T> out_eT;
+  typedef typename T1::pod_type T;
  
-  if(in.aux_out_eT != out_eT(0))
-    {
-    out.init_xform_mt(in.m, priv::functor_cx_scalar_times<T>(in.aux_out_eT));
-    }
-  else
-    {
-    const SpProxy<T1> P(in.m);
-    
-    out.zeros( P.get_n_rows(), P.get_n_cols() );
-    }
+  out.init_xform_mt(in.m, priv::functor_cx_scalar_times<T>(in.aux_out_eT));
  }


@@ -566,4 +547,21 @@ spop_fliplr::apply(SpMat<typename T1::elem_type>& out, const SpOp<T1,spop_fliplr



+template<typename eT, typename T1>
+inline
+void
+spop_replace::apply(SpMat<eT>& out, const mtSpOp<eT, T1, spop_replace>& in)
+  {
+  arma_debug_sigprint();
+  
+  const eT old_val = in.aux;
+  const eT new_val = in.aux_out_eT;
+  
+  out = in.m;
+  
+  out.replace(old_val, new_val);
+  }
+
+
+
 //! @}
@@ -0,0 +1,35 @@
+// SPDX-License-Identifier: Apache-2.0
+// 
+// Copyright 2008-2016 Conrad Sanderson (http://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
+// http://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 spop_omit
+//! @{
+
+
+class spop_omit
+  : public traits_op_col
+  {
+  public:
+  
+  template<typename T1> inline static void apply(SpMat<typename T1::elem_type>& out, const SpOp<T1, spop_omit>& in);
+  
+  template<typename T1, typename functor> inline static void apply_noalias(SpMat<typename T1::elem_type>& out, const SpProxy<T1>& P, functor is_omitted);
+  };
+
+
+//! @}
@@ -0,0 +1,119 @@
+// SPDX-License-Identifier: Apache-2.0
+// 
+// Copyright 2008-2016 Conrad Sanderson (http://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
+// http://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 spop_omit
+//! @{
+
+
+
+template<typename T1>
+inline
+void
+spop_omit::apply(SpMat<typename T1::elem_type>& out, const SpOp<T1, spop_omit>& in)
+  {
+  arma_debug_sigprint();
+  
+  typedef typename T1::elem_type eT;
+  
+  const uword omit_mode = in.aux_uword_a;
+  
+  if(arma_config::fast_math_warn)
+    {
+    if(omit_mode == 1)  { arma_warn(1, "omit_nan(): detection of NaN is not reliable in fast math mode"); }
+    if(omit_mode == 2)  { arma_warn(1, "omit_nonfinite(): detection of non-finite values is not reliable in fast math mode"); }
+    }
+  
+  auto is_omitted_1 = [](const eT& x) -> bool { return arma_isnan(x);       };
+  auto is_omitted_2 = [](const eT& x) -> bool { return arma_isnonfinite(x); };
+  
+  const SpProxy<T1> P(in.m);
+  
+  if(P.is_alias(out))
+    {
+    SpMat<eT> tmp;
+    
+    if(omit_mode == 1)  { spop_omit::apply_noalias(tmp, P, is_omitted_1); }
+    if(omit_mode == 2)  { spop_omit::apply_noalias(tmp, P, is_omitted_2); }
+    
+    out.steal_mem(tmp);
+    }
+  else
+    {
+    if(omit_mode == 1)  { spop_omit::apply_noalias(out, P, is_omitted_1); }
+    if(omit_mode == 2)  { spop_omit::apply_noalias(out, P, is_omitted_2); }
+    }
+  }
+
+
+
+template<typename T1, typename functor>
+inline
+void
+spop_omit::apply_noalias(SpMat<typename T1::elem_type>& out, const SpProxy<T1>& P, functor is_omitted)
+  {
+  arma_debug_sigprint();
+  
+  typedef typename T1::elem_type eT;
+  
+  const uword n_rows        = P.get_n_rows();
+  const uword max_n_nonzero = P.get_n_nonzero();
+  
+  if(max_n_nonzero == 0)  { out.reset(); return; }
+  
+  out.reserve(P.get_n_elem(), 1, max_n_nonzero);
+  
+  typename SpProxy<T1>::const_iterator_type it     = P.begin();
+  typename SpProxy<T1>::const_iterator_type it_end = P.end();
+  
+  uword count = 0;
+  
+  for(; it != it_end; ++it)
+    {
+    const eT val = (*it);
+    
+    if(is_omitted(val) == false)
+      {
+      const uword index = it.row() + it.col()*n_rows;
+      
+      access::rw(out.values[count])      = val;
+      access::rw(out.row_indices[count]) = index;
+      access::rw(out.col_ptrs[1])++;
+      ++count;
+      }
+    }
+  
+  if(count < max_n_nonzero)
+    {
+    if(count <= (max_n_nonzero/2))
+      {
+      out.mem_resize(count);
+      }
+    else
+      {
+      // quick resize without reallocating memory and copying data
+      access::rw(         out.n_nonzero) = count;
+      access::rw(     out.values[count]) = eT(0);
+      access::rw(out.row_indices[count]) = uword(0);
+      }
+    }
+  }
+
+
+
+//! @}
@@ -474,8 +474,8 @@ subview_elem1<eT,T1>::randu()
        eT*   m_mem    = m_local.memptr();
  const uword m_n_elem = m_local.n_elem;
  
-  const unwrap_check_mixed<T1> tmp(a.get_ref(), m_local);
-  const umat& aa = tmp.M;
+  const unwrap_check_mixed<T1> U(a.get_ref(), m_local);
+  const umat& aa = U.M;
  
  if(resolves_to_vector<T1>::no)
    {
@@ -485,28 +485,19 @@ subview_elem1<eT,T1>::randu()
  const uword* aa_mem    = aa.memptr();
  const uword  aa_n_elem = aa.n_elem;
  
-  uword iq,jq;
-  for(iq=0, jq=1; jq < aa_n_elem; iq+=2, jq+=2)
-    {
-    const uword ii = aa_mem[iq];
-    const uword jj = aa_mem[jq];
-    
-    arma_conform_check_bounds( ( (ii >= m_n_elem) || (jj >= m_n_elem) ), "Mat::elem(): index out of bounds" );
-    
-    const eT val1 = eT(arma_rng::randu<eT>());
-    const eT val2 = eT(arma_rng::randu<eT>());
-    
-    m_mem[ii] = val1;
-    m_mem[jj] = val2;
-    }
+  podarray<eT> tmp(aa_n_elem);
  
-  if(iq < aa_n_elem)
+  eT* tmp_mem = tmp.memptr();
+  
+  arma_rng::randu<eT>::fill(tmp_mem, aa_n_elem);
+  
+  for(uword iq=0; iq < aa_n_elem; ++iq)
    {
    const uword ii = aa_mem[iq];
    
-    arma_conform_check_bounds( (ii >= m_n_elem) , "Mat::elem(): index out of bounds" ); 
+    arma_conform_check_bounds( (ii >= m_n_elem), "Mat::elem(): index out of bounds" );
    
-    m_mem[ii] = eT(arma_rng::randu<eT>());
+    m_mem[ii] = tmp_mem[iq];
    }
  }

@@ -524,8 +515,8 @@ subview_elem1<eT,T1>::randn()
        eT*   m_mem    = m_local.memptr();
  const uword m_n_elem = m_local.n_elem;
  
-  const unwrap_check_mixed<T1> tmp(a.get_ref(), m_local);
-  const umat& aa = tmp.M;
+  const unwrap_check_mixed<T1> U(a.get_ref(), m_local);
+  const umat& aa = U.M;
  
  if(resolves_to_vector<T1>::no)
    {
@@ -535,24 +526,19 @@ subview_elem1<eT,T1>::randn()
  const uword* aa_mem    = aa.memptr();
  const uword  aa_n_elem = aa.n_elem;
  
-  uword iq,jq;
-  for(iq=0, jq=1; jq < aa_n_elem; iq+=2, jq+=2)
-    {
-    const uword ii = aa_mem[iq];
-    const uword jj = aa_mem[jq];
-    
-    arma_conform_check_bounds( ( (ii >= m_n_elem) || (jj >= m_n_elem) ), "Mat::elem(): index out of bounds" );
-    
-    arma_rng::randn<eT>::dual_val( m_mem[ii], m_mem[jj] );
-    }
+  podarray<eT> tmp(aa_n_elem);
  
-  if(iq < aa_n_elem)
+  eT* tmp_mem = tmp.memptr();
+  
+  arma_rng::randn<eT>::fill(tmp_mem, aa_n_elem);
+  
+  for(uword iq=0; iq < aa_n_elem; ++iq)
    {
    const uword ii = aa_mem[iq];
    
-    arma_conform_check_bounds( (ii >= m_n_elem) , "Mat::elem(): index out of bounds" ); 
+    arma_conform_check_bounds( (ii >= m_n_elem), "Mat::elem(): index out of bounds" );
    
-    m_mem[ii] = eT(arma_rng::randn<eT>());
+    m_mem[ii] = tmp_mem[iq];
    }
  }

@@ -67,6 +67,8 @@ class subview_elem2 : public Base< eT, subview_elem2<eT,T1,T2> >
  inline void fill(const eT val);
  inline void zeros();
  inline void ones();
+  inline void randu();
+  inline void randn();
  
  inline void operator+= (const eT val);
  inline void operator-= (const eT val);
@@ -299,6 +299,240 @@ subview_elem2<eT,T1,T2>::inplace_op(const Base<eT,expr>& x)



+template<typename eT, typename T1, typename T2>
+inline
+void
+subview_elem2<eT,T1,T2>::randu()
+  {
+  arma_debug_sigprint();
+  
+  Mat<eT>& m_local = const_cast< Mat<eT>& >(m);
+  
+  const uword m_n_rows = m_local.n_rows;
+  const uword m_n_cols = m_local.n_cols;
+  
+  if( (all_rows == false) && (all_cols == false) )
+    {
+    const unwrap_check_mixed<T1> U1(base_ri.get_ref(), m_local);
+    const unwrap_check_mixed<T2> U2(base_ci.get_ref(), m_local);
+    
+    const umat& ri = U1.M;
+    const umat& ci = U2.M;
+    
+    arma_conform_check
+      (
+      ( ((ri.is_vec() == false) && (ri.is_empty() == false)) || ((ci.is_vec() == false) && (ci.is_empty() == false)) ),
+      "Mat::elem(): given object must be a vector"
+      );
+    
+    const uword* ri_mem    = ri.memptr();
+    const uword  ri_n_elem = ri.n_elem;
+    
+    const uword* ci_mem    = ci.memptr();
+    const uword  ci_n_elem = ci.n_elem;
+    
+    podarray<eT> tmp(ri_n_elem);
+    
+    eT* tmp_mem = tmp.memptr();
+    
+    for(uword ci_count=0; ci_count < ci_n_elem; ++ci_count)
+      {
+      const uword col = ci_mem[ci_count];
+      
+      arma_conform_check_bounds( (col >= m_n_cols), "Mat::elem(): index out of bounds" );
+      
+      arma_rng::randu<eT>::fill(tmp_mem, ri_n_elem);
+      
+      for(uword ri_count=0; ri_count < ri_n_elem; ++ri_count)
+        {
+        const uword row = ri_mem[ri_count];
+        
+        arma_conform_check_bounds( (row >= m_n_rows), "Mat::elem(): index out of bounds" );
+        
+        m_local.at(row,col) = tmp_mem[ri_count];
+        }
+      }
+    }
+  else
+  if( (all_rows == true) && (all_cols == false) )
+    {
+    const unwrap_check_mixed<T2> U2(base_ci.get_ref(), m_local);
+    
+    const umat& ci = U2.M;
+    
+    arma_conform_check
+      (
+      ( (ci.is_vec() == false) && (ci.is_empty() == false) ),
+      "Mat::elem(): given object must be a vector"
+      );
+    
+    const uword* ci_mem    = ci.memptr();
+    const uword  ci_n_elem = ci.n_elem;
+    
+    for(uword ci_count=0; ci_count < ci_n_elem; ++ci_count)
+      {
+      const uword col = ci_mem[ci_count];
+      
+      arma_conform_check_bounds( (col >= m_n_cols), "Mat::elem(): index out of bounds" );
+      
+      arma_rng::randu<eT>::fill(m_local.colptr(col), m_n_rows);
+      }
+    }
+  else
+  if( (all_rows == false) && (all_cols == true) )
+    {
+    const unwrap_check_mixed<T1> U1(base_ri.get_ref(), m_local);
+    
+    const umat& ri = U1.M;
+    
+    arma_conform_check
+      (
+      ( (ri.is_vec() == false) && (ri.is_empty() == false) ),
+      "Mat::elem(): given object must be a vector"
+      );
+    
+    const uword* ri_mem    = ri.memptr();
+    const uword  ri_n_elem = ri.n_elem;
+    
+    podarray<eT> tmp(ri_n_elem);
+    
+    eT* tmp_mem = tmp.memptr();
+    
+    for(uword col=0; col < m_n_cols; ++col)
+      {
+      arma_rng::randu<eT>::fill(tmp_mem, ri_n_elem);
+      
+      for(uword ri_count=0; ri_count < ri_n_elem; ++ri_count)
+        {
+        const uword row = ri_mem[ri_count];
+        
+        arma_conform_check_bounds( (row >= m_n_rows), "Mat::elem(): index out of bounds" );
+      
+        m_local.at(row,col) = tmp_mem[ri_count];
+        }
+      }
+    }
+  }
+
+
+
+template<typename eT, typename T1, typename T2>
+inline
+void
+subview_elem2<eT,T1,T2>::randn()
+  {
+  arma_debug_sigprint();
+  
+  Mat<eT>& m_local = const_cast< Mat<eT>& >(m);
+  
+  const uword m_n_rows = m_local.n_rows;
+  const uword m_n_cols = m_local.n_cols;
+  
+  if( (all_rows == false) && (all_cols == false) )
+    {
+    const unwrap_check_mixed<T1> U1(base_ri.get_ref(), m_local);
+    const unwrap_check_mixed<T2> U2(base_ci.get_ref(), m_local);
+    
+    const umat& ri = U1.M;
+    const umat& ci = U2.M;
+    
+    arma_conform_check
+      (
+      ( ((ri.is_vec() == false) && (ri.is_empty() == false)) || ((ci.is_vec() == false) && (ci.is_empty() == false)) ),
+      "Mat::elem(): given object must be a vector"
+      );
+    
+    const uword* ri_mem    = ri.memptr();
+    const uword  ri_n_elem = ri.n_elem;
+    
+    const uword* ci_mem    = ci.memptr();
+    const uword  ci_n_elem = ci.n_elem;
+    
+    podarray<eT> tmp(ri_n_elem);
+    
+    eT* tmp_mem = tmp.memptr();
+    
+    for(uword ci_count=0; ci_count < ci_n_elem; ++ci_count)
+      {
+      const uword col = ci_mem[ci_count];
+      
+      arma_conform_check_bounds( (col >= m_n_cols), "Mat::elem(): index out of bounds" );
+      
+      arma_rng::randn<eT>::fill(tmp_mem, ri_n_elem);
+      
+      for(uword ri_count=0; ri_count < ri_n_elem; ++ri_count)
+        {
+        const uword row = ri_mem[ri_count];
+        
+        arma_conform_check_bounds( (row >= m_n_rows), "Mat::elem(): index out of bounds" );
+        
+        m_local.at(row,col) = tmp_mem[ri_count];
+        }
+      }
+    }
+  else
+  if( (all_rows == true) && (all_cols == false) )
+    {
+    const unwrap_check_mixed<T2> U2(base_ci.get_ref(), m_local);
+    
+    const umat& ci = U2.M;
+    
+    arma_conform_check
+      (
+      ( (ci.is_vec() == false) && (ci.is_empty() == false) ),
+      "Mat::elem(): given object must be a vector"
+      );
+    
+    const uword* ci_mem    = ci.memptr();
+    const uword  ci_n_elem = ci.n_elem;
+    
+    for(uword ci_count=0; ci_count < ci_n_elem; ++ci_count)
+      {
+      const uword col = ci_mem[ci_count];
+      
+      arma_conform_check_bounds( (col >= m_n_cols), "Mat::elem(): index out of bounds" );
+      
+      arma_rng::randn<eT>::fill(m_local.colptr(col), m_n_rows);
+      }
+    }
+  else
+  if( (all_rows == false) && (all_cols == true) )
+    {
+    const unwrap_check_mixed<T1> U1(base_ri.get_ref(), m_local);
+    
+    const umat& ri = U1.M;
+    
+    arma_conform_check
+      (
+      ( (ri.is_vec() == false) && (ri.is_empty() == false) ),
+      "Mat::elem(): given object must be a vector"
+      );
+    
+    const uword* ri_mem    = ri.memptr();
+    const uword  ri_n_elem = ri.n_elem;
+    
+    podarray<eT> tmp(ri_n_elem);
+    
+    eT* tmp_mem = tmp.memptr();
+    
+    for(uword col=0; col < m_n_cols; ++col)
+      {
+      arma_rng::randn<eT>::fill(tmp_mem, ri_n_elem);
+      
+      for(uword ri_count=0; ri_count < ri_n_elem; ++ri_count)
+        {
+        const uword row = ri_mem[ri_count];
+        
+        arma_conform_check_bounds( (row >= m_n_rows), "Mat::elem(): index out of bounds" );
+      
+        m_local.at(row,col) = tmp_mem[ri_count];
+        }
+      }
+    }
+  }
+
+
+
 //
 //

@@ -58,8 +58,8 @@ guess_sympd_worker(const Mat<eT>& A)
    {
    const eT A_jj = A_col[j];
    
-    if(              A_jj  <= eT(0))  { return false; }
-    if(arma_isfinite(A_jj) == false)  { return false; }
+    if(        A_jj  <= eT(0))  { return false; }
+    if(arma_isnonfinite(A_jj))  { return false; }
    
    if(A_jj >= tol)  { diag_below_tol = false; }
    
@@ -147,8 +147,8 @@ guess_sympd_worker(const Mat<eT>& A)
    const  T  A_jj_rabs =  std::abs(A_jj_r);
    const  T  A_jj_iabs =  std::abs(A_jj_i);
    
-    if(              A_jj_r  <= T(0) )  { return false; }  // real should be positive
-    if(arma_isfinite(A_jj_r) == false)  { return false; }
+    if(        A_jj_r  <= T(0) )  { return false; }  // real should be positive
+    if(arma_isnonfinite(A_jj_r))  { return false; }
    
    if(A_jj_iabs > tol      )  { return false; }  // imag should be approx zero
    if(A_jj_iabs > A_jj_rabs)  { return false; }  // corner case: real and imag are close to zero, and imag is dominant
@@ -164,7 +164,7 @@ guess_sympd_worker(const Mat<eT>& A)
  
  const T square_max_diag = max_diag * max_diag;
  
-  if(arma_isfinite(square_max_diag) == false)  { return false; }
+  if(arma_isnonfinite(square_max_diag))  { return false; }
  
  A_col = A_mem;
  
@@ -188,7 +188,7 @@ guess_sympd_worker(const Mat<eT>& A)
      // avoid using std::abs(), as that is time consuming due to division and std::sqrt()
      const T square_A_ij_abs = (A_ij_real * A_ij_real) + (A_ij_imag * A_ij_imag);
      
-      if(arma_isfinite(square_A_ij_abs) == false)  { return false; }
+      if(arma_isnonfinite(square_A_ij_abs))  { return false; }
      
      if(square_A_ij_abs >= square_max_diag)  { return false; }
      
@@ -285,7 +285,7 @@ is_approx_sym_worker(const Mat<eT>& A)
    {
    const eT A_jj = A_col[j];
    
-    if(arma_isfinite(A_jj) == false)  { return false; }
+    if(arma_isnonfinite(A_jj))  { return false; }
    
    if(std::abs(A_jj) >= tol)  { diag_below_tol = false; }
    
@@ -359,7 +359,7 @@ is_approx_sym_worker(const Mat<eT>& A)
    if(A_jj_iabs > tol      )  { return false; }  // imag should be approx zero
    if(A_jj_iabs > A_jj_rabs)  { return false; }  // corner case: real and imag are close to zero, and imag is dominant
    
-    if(arma_isfinite(A_jj_r) == false)  { return false; }
+    if(arma_isnonfinite(A_jj_r))  { return false; }
    
    if(A_jj_rabs >= tol)  { diag_below_tol = false; }
    
@@ -1485,6 +1485,28 @@ namespace lapack
    #endif
    }
  
+  
+  
+  template<typename eT>
+  inline
+  void
+  gebal(const char* job, const blas_int* n, eT* a, const blas_int* lda, blas_int* ilo, blas_int* ihi, typename get_pod_type<eT>::result* scale, blas_int* info)
+    {
+    arma_type_check(( is_supported_blas_type<eT>::value == false ));
+    
+    #if defined(ARMA_USE_FORTRAN_HIDDEN_ARGS)
+           if(    is_float<eT>::value)  { typedef float  pod_T; typedef float    T; arma_fortran(arma_sgebal)(job, n, (T*)a, lda, ilo, ihi, (pod_T*)scale, info, 1); }
+      else if(   is_double<eT>::value)  { typedef double pod_T; typedef double   T; arma_fortran(arma_dgebal)(job, n, (T*)a, lda, ilo, ihi, (pod_T*)scale, info, 1); }
+      else if( is_cx_float<eT>::value)  { typedef float  pod_T; typedef blas_cxf T; arma_fortran(arma_cgebal)(job, n, (T*)a, lda, ilo, ihi, (pod_T*)scale, info, 1); }
+      else if(is_cx_double<eT>::value)  { typedef double pod_T; typedef blas_cxd T; arma_fortran(arma_zgebal)(job, n, (T*)a, lda, ilo, ihi, (pod_T*)scale, info, 1); }
+    #else
+           if(    is_float<eT>::value)  { typedef float  pod_T; typedef float    T; arma_fortran(arma_sgebal)(job, n, (T*)a, lda, ilo, ihi, (pod_T*)scale, info); }
+      else if(   is_double<eT>::value)  { typedef double pod_T; typedef double   T; arma_fortran(arma_dgebal)(job, n, (T*)a, lda, ilo, ihi, (pod_T*)scale, info); }
+      else if( is_cx_float<eT>::value)  { typedef float  pod_T; typedef blas_cxf T; arma_fortran(arma_cgebal)(job, n, (T*)a, lda, ilo, ihi, (pod_T*)scale, info); }
+      else if(is_cx_double<eT>::value)  { typedef double pod_T; typedef blas_cxd T; arma_fortran(arma_zgebal)(job, n, (T*)a, lda, ilo, ihi, (pod_T*)scale, info); }
+    #endif
+    }
+  
  }


@@ -26,7 +26,7 @@ namespace superlu
  template<typename eT>
  inline
  void
-  gssv(superlu_options_t* options, SuperMatrix* A, int* perm_c, int* perm_r, SuperMatrix* L, SuperMatrix* U, SuperMatrix* B, SuperLUStat_t* stat, int* info)
+  gssv(superlu_options_t* options, SuperMatrix* A, int* perm_c, int* perm_r, SuperMatrix* L, SuperMatrix* U, SuperMatrix* B, SuperLUStat_t* stat, superlu::int_t* info)
    {
    arma_type_check(( is_supported_blas_type<eT>::value == false ));
    
@@ -63,11 +63,11 @@ namespace superlu
        int* etree, char* equed,
        typename get_pod_type<eT>::result* R, typename get_pod_type<eT>::result* C,
        SuperMatrix* L, SuperMatrix* U,
-        void* work, int lwork,
+        void* work, superlu::int_t lwork,
        SuperMatrix* B, SuperMatrix* X,
        typename get_pod_type<eT>::result* rpg, typename get_pod_type<eT>::result* rcond,
        typename get_pod_type<eT>::result* ferr, typename get_pod_type<eT>::result* berr,
-        GlobalLU_t* glu, mem_usage_t* mu, SuperLUStat_t* stat, int* info
+        GlobalLU_t* glu, mem_usage_t* mu, SuperLUStat_t* stat, superlu::int_t* info
       )
    {
    arma_type_check(( is_supported_blas_type<eT>::value == false ));
@@ -105,11 +105,11 @@ namespace superlu
  gstrf(superlu_options_t* options,
        SuperMatrix* A,
        int relax,
-        int panel_size, int *etree,
-        void  *work,  int  lwork,
+        int panel_size, int* etree,
+        void* work, superlu::int_t lwork,
        int* perm_c, int* perm_r,
        SuperMatrix* L, SuperMatrix* U,
-        GlobalLU_t* Glu, SuperLUStat_t* stat, int* info
+        GlobalLU_t* Glu, SuperLUStat_t* stat, superlu::int_t* info
       )
    {
    arma_type_check(( is_supported_blas_type<eT>::value == false ));
@@ -15,10 +15,10 @@
 </li>
 </ul>
 <ul>
-<li>Article with overview of functionality in Armadillo
+<li>Paper with overview of functionality in Armadillo
 <ul>
-<li><a href="https://arma.sourceforge.net/armadillo_joss_2016.pdf">armadillo_joss_2016.pdf</a> &nbsp; (online copy)</li>
-<li><a href="armadillo_joss_2016.pdf">armadillo_joss_2016.pdf</a> &nbsp; (local copy)</li>
+<li><a href="https://arma.sourceforge.net/armadillo_iccae_2025.pdf">armadillo_iccae_2025.pdf</a> &nbsp; (online copy)</li>
+<li><a href="armadillo_iccae_2025.pdf">armadillo_iccae_2025.pdf</a> &nbsp; (local copy)</li>
 </ul>
 </li>
 </ul>
@@ -27,10 +27,6 @@
 <ul>
 <li><a href="https://arma.sourceforge.net/armadillo_mca_2019.pdf">armadillo_mca_2019.pdf</a> &nbsp; (online copy)</li>
 <li><a href="armadillo_mca_2019.pdf">armadillo_mca_2019.pdf</a> &nbsp; (local copy)</li>
-<!--
-<li><a href="https://arma.sourceforge.net/armadillo_lncs_2018.pdf">armadillo_lncs_2018.pdf</a> &nbsp; (online copy)</li>
-<li><a href="armadillo_lncs_2018.pdf">armadillo_lncs_2018.pdf</a> &nbsp; (local copy)</li>
-->
 </ul>
 </li>
 </ul>
@@ -51,14 +47,6 @@
 </li>
 </ul>
 <ul>
-<li>Overview technical report (old)
-<ul>
-<li><a href="https://arma.sourceforge.net/armadillo_nicta_2010.pdf">armadillo_nicta_2010.pdf</a> &nbsp; (online copy)</li>
-<li><a href="armadillo_nicta_2010.pdf">armadillo_nicta_2010.pdf</a> &nbsp; (local copy)</li>
-</ul>
-</li>
-</ul>
-<ul>
 <li>Article on integration of Armadillo with R:
 <ul>
 <li><a href="https://arma.sourceforge.net/armadillo_rcpp_2014.pdf">armadillo_rcpp_2014.pdf</a> &nbsp; (online copy)</li>
@@ -1383,6 +1383,28 @@ extern "C"
      arma_fortran_sans_prefix(arma_zhecon)(uplo, n, a, lda, ipiv, anorm, rcond, work, info);
      }
    
+    
+    
+    void arma_fortran_with_prefix(arma_sgebal)(const char* job, const blas_int* n,    float* a, const blas_int* lda, blas_int* ilo, blas_int* ihi,  float* scale, blas_int* info)
+      {
+      arma_fortran_sans_prefix(arma_sgebal)(job, n, a, lda, ilo, ihi, scale, info);
+      }
+    
+    void arma_fortran_with_prefix(arma_dgebal)(const char* job, const blas_int* n,   double* a, const blas_int* lda, blas_int* ilo, blas_int* ihi, double* scale, blas_int* info)
+      {
+      arma_fortran_sans_prefix(arma_dgebal)(job, n, a, lda, ilo, ihi, scale, info);
+      }
+    
+    void arma_fortran_with_prefix(arma_cgebal)(const char* job, const blas_int* n, blas_cxf* a, const blas_int* lda, blas_int* ilo, blas_int* ihi,  float* scale, blas_int* info)
+      {
+      arma_fortran_sans_prefix(arma_cgebal)(job, n, a, lda, ilo, ihi, scale, info);
+      }
+    
+    void arma_fortran_with_prefix(arma_zgebal)(const char* job, const blas_int* n, blas_cxd* a, const blas_int* lda, blas_int* ilo, blas_int* ihi, double* scale, blas_int* info)
+      {
+      arma_fortran_sans_prefix(arma_zgebal)(job, n, a, lda, ilo, ihi, scale, info);
+      }
+    
  #endif
  
  
@@ -1458,22 +1480,22 @@ extern "C"
  
  #if defined(ARMA_USE_SUPERLU)
    
-    void wrapper_sgssv(superlu::superlu_options_t* a, superlu::SuperMatrix* b, int* c, int* d, superlu::SuperMatrix* e, superlu::SuperMatrix* f, superlu::SuperMatrix* g, superlu::SuperLUStat_t* h, int* i)
+    void wrapper_sgssv(superlu::superlu_options_t* a, superlu::SuperMatrix* b, int* c, int* d, superlu::SuperMatrix* e, superlu::SuperMatrix* f, superlu::SuperMatrix* g, superlu::SuperLUStat_t* h, superlu::int_t* i)
      {
      sgssv(a,b,c,d,e,f,g,h,i);
      }
    
-    void wrapper_dgssv(superlu::superlu_options_t* a, superlu::SuperMatrix* b, int* c, int* d, superlu::SuperMatrix* e, superlu::SuperMatrix* f, superlu::SuperMatrix* g, superlu::SuperLUStat_t* h, int* i)
+    void wrapper_dgssv(superlu::superlu_options_t* a, superlu::SuperMatrix* b, int* c, int* d, superlu::SuperMatrix* e, superlu::SuperMatrix* f, superlu::SuperMatrix* g, superlu::SuperLUStat_t* h, superlu::int_t* i)
      {
      dgssv(a,b,c,d,e,f,g,h,i);
      }
    
-    void wrapper_cgssv(superlu::superlu_options_t* a, superlu::SuperMatrix* b, int* c, int* d, superlu::SuperMatrix* e, superlu::SuperMatrix* f, superlu::SuperMatrix* g, superlu::SuperLUStat_t* h, int* i)
+    void wrapper_cgssv(superlu::superlu_options_t* a, superlu::SuperMatrix* b, int* c, int* d, superlu::SuperMatrix* e, superlu::SuperMatrix* f, superlu::SuperMatrix* g, superlu::SuperLUStat_t* h, superlu::int_t* i)
      {
      cgssv(a,b,c,d,e,f,g,h,i);
      }
    
-    void wrapper_zgssv(superlu::superlu_options_t* a, superlu::SuperMatrix* b, int* c, int* d, superlu::SuperMatrix* e, superlu::SuperMatrix* f, superlu::SuperMatrix* g, superlu::SuperLUStat_t* h, int* i)
+    void wrapper_zgssv(superlu::superlu_options_t* a, superlu::SuperMatrix* b, int* c, int* d, superlu::SuperMatrix* e, superlu::SuperMatrix* f, superlu::SuperMatrix* g, superlu::SuperLUStat_t* h, superlu::int_t* i)
      {
      zgssv(a,b,c,d,e,f,g,h,i);
      }
@@ -1481,22 +1503,22 @@ extern "C"
    
    
    
-    void wrapper_sgssvx(superlu::superlu_options_t* a, superlu::SuperMatrix* b, int* c, int* d, int* e, char* f,  float* g,  float* h, superlu::SuperMatrix* i, superlu::SuperMatrix* j, void* k, int l, superlu::SuperMatrix* m, superlu::SuperMatrix* n,  float* o,  float* p,  float* q,  float* r, superlu::GlobalLU_t* s, superlu::mem_usage_t* t, superlu::SuperLUStat_t* u, int* v)
+    void wrapper_sgssvx(superlu::superlu_options_t* a, superlu::SuperMatrix* b, int* c, int* d, int* e, char* f,  float* g,  float* h, superlu::SuperMatrix* i, superlu::SuperMatrix* j, void* k, superlu::int_t l, superlu::SuperMatrix* m, superlu::SuperMatrix* n,  float* o,  float* p,  float* q,  float* r, superlu::GlobalLU_t* s, superlu::mem_usage_t* t, superlu::SuperLUStat_t* u, superlu::int_t* v)
      {
      sgssvx(a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p, q, r, s, t, u, v);
      }
    
-    void wrapper_dgssvx(superlu::superlu_options_t* a, superlu::SuperMatrix* b, int* c, int* d, int* e, char* f, double* g, double* h, superlu::SuperMatrix* i, superlu::SuperMatrix* j, void* k, int l, superlu::SuperMatrix* m, superlu::SuperMatrix* n, double* o, double* p, double* q, double* r, superlu::GlobalLU_t* s, superlu::mem_usage_t* t, superlu::SuperLUStat_t* u, int* v)
+    void wrapper_dgssvx(superlu::superlu_options_t* a, superlu::SuperMatrix* b, int* c, int* d, int* e, char* f, double* g, double* h, superlu::SuperMatrix* i, superlu::SuperMatrix* j, void* k, superlu::int_t l, superlu::SuperMatrix* m, superlu::SuperMatrix* n, double* o, double* p, double* q, double* r, superlu::GlobalLU_t* s, superlu::mem_usage_t* t, superlu::SuperLUStat_t* u, superlu::int_t* v)
      {
      dgssvx(a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p, q, r, s, t, u, v);
      }
    
-    void wrapper_cgssvx(superlu::superlu_options_t* a, superlu::SuperMatrix* b, int* c, int* d, int* e, char* f,  float* g,  float* h, superlu::SuperMatrix* i, superlu::SuperMatrix* j, void* k, int l, superlu::SuperMatrix* m, superlu::SuperMatrix* n,  float* o,  float* p,  float* q,  float* r, superlu::GlobalLU_t* s, superlu::mem_usage_t* t, superlu::SuperLUStat_t* u, int* v)
+    void wrapper_cgssvx(superlu::superlu_options_t* a, superlu::SuperMatrix* b, int* c, int* d, int* e, char* f,  float* g,  float* h, superlu::SuperMatrix* i, superlu::SuperMatrix* j, void* k, superlu::int_t l, superlu::SuperMatrix* m, superlu::SuperMatrix* n,  float* o,  float* p,  float* q,  float* r, superlu::GlobalLU_t* s, superlu::mem_usage_t* t, superlu::SuperLUStat_t* u, superlu::int_t* v)
      {
      cgssvx(a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p, q, r, s, t, u, v);
      }
    
-    void wrapper_zgssvx(superlu::superlu_options_t* a, superlu::SuperMatrix* b, int* c, int* d, int* e, char* f, double* g, double* h, superlu::SuperMatrix* i, superlu::SuperMatrix* j, void* k, int l, superlu::SuperMatrix* m, superlu::SuperMatrix* n, double* o, double* p, double* q, double* r, superlu::GlobalLU_t* s, superlu::mem_usage_t* t, superlu::SuperLUStat_t* u, int* v)
+    void wrapper_zgssvx(superlu::superlu_options_t* a, superlu::SuperMatrix* b, int* c, int* d, int* e, char* f, double* g, double* h, superlu::SuperMatrix* i, superlu::SuperMatrix* j, void* k, superlu::int_t l, superlu::SuperMatrix* m, superlu::SuperMatrix* n, double* o, double* p, double* q, double* r, superlu::GlobalLU_t* s, superlu::mem_usage_t* t, superlu::SuperLUStat_t* u, superlu::int_t* v)
      {
      zgssvx(a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p, q, r, s, t, u, v);
      }
@@ -1504,22 +1526,22 @@ extern "C"
    
    
    
-    void wrapper_sgstrf(superlu::superlu_options_t* a, superlu::SuperMatrix* b, int d, int e, int* f, void* g, int h, int* i, int* j, superlu::SuperMatrix* k, superlu::SuperMatrix* l, superlu::GlobalLU_t* m, superlu::SuperLUStat_t* n, int* o)
+    void wrapper_sgstrf(superlu::superlu_options_t* a, superlu::SuperMatrix* b, int d, int e, int* f, void* g, superlu::int_t h, int* i, int* j, superlu::SuperMatrix* k, superlu::SuperMatrix* l, superlu::GlobalLU_t* m, superlu::SuperLUStat_t* n, superlu::int_t* o)
      {
      sgstrf(a, b, d, e, f, g, h, i, j, k, l, m, n, o);
      }
    
-    void wrapper_dgstrf(superlu::superlu_options_t* a, superlu::SuperMatrix* b, int d, int e, int* f, void* g, int h, int* i, int* j, superlu::SuperMatrix* k, superlu::SuperMatrix* l, superlu::GlobalLU_t* m, superlu::SuperLUStat_t* n, int* o)
+    void wrapper_dgstrf(superlu::superlu_options_t* a, superlu::SuperMatrix* b, int d, int e, int* f, void* g, superlu::int_t h, int* i, int* j, superlu::SuperMatrix* k, superlu::SuperMatrix* l, superlu::GlobalLU_t* m, superlu::SuperLUStat_t* n, superlu::int_t* o)
      {
      dgstrf(a, b, d, e, f, g, h, i, j, k, l, m, n, o);
      }
    
-    void wrapper_cgstrf(superlu::superlu_options_t* a, superlu::SuperMatrix* b, int d, int e, int* f, void* g, int h, int* i, int* j, superlu::SuperMatrix* k, superlu::SuperMatrix* l, superlu::GlobalLU_t* m, superlu::SuperLUStat_t* n, int* o)
+    void wrapper_cgstrf(superlu::superlu_options_t* a, superlu::SuperMatrix* b, int d, int e, int* f, void* g, superlu::int_t h, int* i, int* j, superlu::SuperMatrix* k, superlu::SuperMatrix* l, superlu::GlobalLU_t* m, superlu::SuperLUStat_t* n, superlu::int_t* o)
      {
      cgstrf(a, b, d, e, f, g, h, i, j, k, l, m, n, o);
      }
    
-    void wrapper_zgstrf(superlu::superlu_options_t* a, superlu::SuperMatrix* b, int d, int e, int* f, void* g, int h, int* i, int* j, superlu::SuperMatrix* k, superlu::SuperMatrix* l, superlu::GlobalLU_t* m, superlu::SuperLUStat_t* n, int* o)
+    void wrapper_zgstrf(superlu::superlu_options_t* a, superlu::SuperMatrix* b, int d, int e, int* f, void* g, superlu::int_t h, int* i, int* j, superlu::SuperMatrix* k, superlu::SuperMatrix* l, superlu::GlobalLU_t* m, superlu::SuperLUStat_t* n, superlu::int_t* o)
      {
      zgstrf(a, b, d, e, f, g, h, i, j, k, l, m, n, o);
      }
--- a/Show More
+++ b/Show More
Author	SHA1	Message	Date
conrad	98c13f4bc7	patch bump	2025-08-08 13:56:43 +10:00
conrad	92cd541f57	move conformance check to avoid doing two checks	2025-08-05 23:51:07 +10:00
conrad	f002bffd96	simplified fallback in handling of op_sum specialisations	2025-08-05 23:43:29 +10:00
conrad	cba9c87169	fix	2025-08-05 12:39:40 +10:00
conrad	12628de147	patch bump	2025-07-21 11:51:08 +10:00
conrad	16d102510f	show sylvester() instead of deprecatated syl()	2025-07-21 11:40:40 +10:00
conrad	87626c836a	document ARMA_DONT_OPTIMISE_INVEXPR and ARMA_DONT_OPTIMISE_POWEXPR	2025-07-19 01:30:05 +10:00
conrad	6e7f5c3601	add macro to disable optimisation of pow()	2025-07-16 13:21:36 +10:00
conrad	7ac92e7999	restore compile time deprecation message	2025-07-14 14:50:51 +10:00
conrad	fc11d3683d	restore compile time deprecation message	2025-07-14 14:50:03 +10:00
conrad	3eb25b1b16	allow to optionally disable optimised handling of pow()	2025-07-14 14:36:22 +10:00
conrad	a197882af1	better wording	2025-07-12 18:49:12 +10:00
conrad	0dd9ac9a02	better wording	2025-07-12 18:48:09 +10:00
conradsnicta	5d62ce315c	Merge branch 'seb/add_an_ignore' into '14.6.x' Adds the missing argument Xst to the 'arma_ignore list' See merge request conradsnicta/armadillo-code!142	2025-07-12 01:15:44 +00:00
Seb James	51d56e7c59	Adds the missing argument Xst to the 'arma_ignore list'	2025-07-12 01:15:44 +00:00
conrad	17db28c133	remove outdated comments	2025-07-07 21:46:12 +10:00
conrad	b801feb467	simplifications	2025-07-07 20:29:59 +10:00
conrad	ec95397bc4	cleanup	2025-07-07 15:45:57 +10:00
conrad	851828d427	speedup	2025-07-05 13:03:45 +10:00
conrad	1192e1ac19	better wording	2025-07-04 13:31:21 +10:00
conrad	2293bf45cc	cleanup warnings	2025-07-04 13:24:56 +10:00
conrad	bb7ed84ded	expand deprecation message with refactoring suggestion	2025-07-03 15:48:32 +10:00
conrad	d1305141d2	version bump	2025-07-02 16:46:17 +10:00
conrad	f3ed7ace06	mark as RC2	2025-07-01 21:48:15 +10:00
conrad	ed8d293f93	fix	2025-07-01 12:35:50 +10:00
conrad	ce007b7164	restore previous setup	2025-06-30 19:23:14 +10:00
conrad	b2664f138b	restore previous setup	2025-06-30 19:22:21 +10:00
conrad	300bd731ca	cleanup	2025-06-30 19:20:10 +10:00
conrad	9843393a7a	cleanup	2025-06-30 19:18:03 +10:00
conrad	98a885b9c2	mark as RC1	2025-06-30 19:16:06 +10:00
conrad	0b4542fc88	restore previous setup	2025-06-30 19:15:02 +10:00
conrad	5765ef2909	cleanup	2025-06-30 12:31:26 +10:00
conrad	fd2b724040	rename accumulator variable for clarity	2025-06-29 23:06:46 +10:00
conrad	e7bd7e4e18	return early if there are no non-zero elements	2025-06-29 22:52:33 +10:00
conrad	a902935f11	avoid going through SpProxy unless necessary	2025-06-29 00:38:05 +10:00
conrad	1907b8799e	fix debug messages	2025-06-29 00:28:42 +10:00
conrad	9782878612	minor optimisation	2025-06-29 00:22:24 +10:00
conrad	b94080a168	specialised handling of spop_omit	2025-06-29 00:10:36 +10:00
conrad	ada1e18e45	split out handling of spop_square into separate function	2025-06-28 19:16:24 +10:00
conrad	bd26d90563	fix example	2025-06-28 15:40:57 +10:00
conrad	447d66758e	simplification	2025-06-28 14:21:09 +10:00
conrad	d242d1d3d7	add internal link	2025-06-28 14:13:26 +10:00
conrad	e7e189c012	add internal links	2025-06-28 13:46:08 +10:00
conrad	256e1d6e88	use consistent wording for infinity	2025-06-28 13:15:49 +10:00
conrad	6952ca11d5	add caveat	2025-06-28 13:09:27 +10:00
conrad	4ec76d5211	simplifications	2025-06-28 13:03:54 +10:00
conrad	05670bfc64	better wording	2025-06-28 12:56:43 +10:00
conrad	29f4f10118	emit warnings when using fast math mode	2025-06-24 18:48:44 +10:00
conrad	e5663df8fa	specialised handling for op_omit_cube	2025-06-24 18:42:01 +10:00
conrad	cefb070832	initial implementation of op_omit_cube	2025-06-24 18:24:26 +10:00
conrad	2f9829c777	better wording	2025-06-24 12:00:48 +10:00
conrad	78dbddf2e9	emit warnings when using fast math mode	2025-06-24 11:05:36 +10:00
conrad	9318a5df61	ensure consistent expression style	2025-06-24 10:55:54 +10:00
conrad	8e65bd7d21	fix wording	2025-06-24 01:02:35 +10:00
conrad	42022c75e3	initial docs for omit_nan() and omit_nonfinite()	2025-06-23 14:18:33 +10:00
conrad	869c46e39b	remove elem_opts	2025-06-23 12:59:12 +10:00
conrad	3199376bbf	simplify 'arma_isfinite(x) == false' to 'arma_isnonfinite(x)'	2025-06-20 13:03:29 +10:00
conrad	a93c67cc34	update deprecation comments	2025-06-19 16:55:03 +10:00
conrad	753e1e23c7	initial implementation of spop_omit	2025-06-19 14:39:25 +10:00
conrad	021403a5c1	specialised handling of op_omit	2025-06-19 13:25:50 +10:00
conrad	402a453065	initial implementation of omit_nan() and omit_nonfinite()	2025-06-17 21:30:58 +10:00
conrad	5b05fb843d	initial implementation of omit_nan() and omit_nonfinite()	2025-06-17 21:30:38 +10:00
conrad	772addf5c1	restore previous behaviour	2025-06-12 19:04:38 +10:00
conrad	c08e594437	fix	2025-06-12 14:21:14 +10:00
conrad	9797ad1cad	simplifications	2025-06-12 14:18:08 +10:00
conrad	b260cdaa58	remove detection of NaN from plain versions of median()	2025-06-12 14:12:32 +10:00
conrad	e83cb29013	add median() variant to handle elem_opts	2025-06-12 13:55:48 +10:00
conrad	48e19f8e50	fix	2025-06-11 21:45:30 +10:00
conrad	4e02904b18	add stddev() variant to handle elem_opts	2025-06-11 19:44:43 +10:00
conrad	bd1bd14831	fix	2025-06-11 19:21:52 +10:00
conrad	74feb70896	simplifications	2025-06-11 19:08:07 +10:00
conrad	48a33a1fc3	use original option name for now	2025-06-11 15:35:14 +10:00
conrad	36a5188083	minor cleanup	2025-06-11 15:22:03 +10:00
conrad	0731de2367	simplifications	2025-06-11 15:04:22 +10:00
conrad	91b5f68bed	fix	2025-06-11 13:45:41 +10:00
conrad	a290735187	simplifications	2025-06-11 13:32:47 +10:00
conrad	4b4ab0ced8	add var() variant to handle elem_opts	2025-06-11 13:27:22 +10:00
conrad	7faec6fdce	simplifications	2025-06-10 19:30:19 +10:00
conrad	696fd719e8	simplify use of podarray	2025-06-10 19:23:56 +10:00
conrad	aa0c021737	simplifications	2025-06-10 19:00:27 +10:00
conrad	682790d0bc	started on op_var_omit	2025-06-10 11:39:26 +10:00
conrad	47a0e466e7	add TODO	2025-06-10 11:29:23 +10:00
conrad	8842cdc739	directly handle corner cases	2025-06-09 22:14:08 +10:00
conrad	870414aa06	simplifications	2025-06-09 21:48:48 +10:00
conrad	7b03681ce7	simplifications	2025-06-09 20:10:50 +10:00
conrad	a80f63501b	expand cube versions of accu() and sum() to handle elem_opts	2025-06-08 21:50:29 +10:00
conrad	6556e4ae04	rearrange	2025-06-08 20:45:28 +10:00
conrad	648527eb1b	rearrange	2025-06-08 19:27:36 +10:00
conrad	4a848124fa	expand sparse versions of accu() and sum() to handle elem_opts	2025-06-08 15:24:39 +10:00
conrad	39d01b483e	remove unused typedef	2025-06-08 15:21:03 +10:00
conrad	03e92c43ba	expand op_mean_omit to general matrices	2025-06-08 01:23:11 +10:00
conrad	bdf448e8f0	fix	2025-06-07 21:20:39 +10:00
conrad	d0dfac1a32	simplification	2025-06-07 20:29:29 +10:00
conrad	d057c221d8	simplifications + fix handling of inf corner case	2025-06-07 18:10:33 +10:00
conrad	83bb1e7937	add TODOs	2025-06-06 13:23:48 +10:00
conrad	bdc82c2c76	initial work on op_mean_omit	2025-06-06 13:09:39 +10:00
conrad	3cfec929f6	use arma_isnonfinite() for semantic simplification	2025-06-06 13:04:45 +10:00
conrad	74720f366d	expand sum() with omit_mode	2025-06-05 13:55:41 +10:00
conrad	65012ac724	use more flexible omit_mode	2025-06-05 12:26:31 +10:00
conrad	7b0ce7088e	add accu() variant that can omit nan values	2025-06-04 13:18:03 +10:00
conrad	d2930a2763	fix	2025-06-04 13:01:05 +10:00
conrad	9ef5f44c53	expand deprecation message with refactoring suggestion	2025-06-03 14:16:25 +10:00
conrad	6b79e17ec3	use batch random generation to avoid locking rng for each random number	2025-06-02 17:02:33 +10:00
conrad	eaf7825044	add randu() and randn()	2025-06-02 16:55:59 +10:00
conrad	365eca0c16	use batch random generation to avoid locking rng for each random number	2025-06-02 16:24:06 +10:00
conrad	096591bafb	add TODOs	2025-06-02 13:58:56 +10:00
conrad	68e9cf2a01	use batch random generation to avoid locking rng for each random number	2025-06-02 13:34:10 +10:00
conrad	36a2e4eed4	use batch random generation to avoid locking rng for each random number	2025-06-02 13:19:08 +10:00
conrad	7e7fff181b	always use stable sort for sorting indices	2025-05-29 20:01:54 +10:00
conrad	91319547c4	use consistent layout	2025-05-29 19:44:12 +10:00
conrad	6288784171	resolve complex number ties	2025-05-29 19:40:55 +10:00
conrad	74eb504435	use stable sort for sort_index()	2025-05-29 18:17:43 +10:00
conrad	7e63c9f48c	add docs for standalone replace()	2025-05-28 13:59:28 +10:00
conrad	5560281526	add standalone replace()	2025-05-24 16:23:03 +10:00
conrad	9d04feb8c3	add dual aux constructor	2025-05-24 16:21:17 +10:00
conrad	47d576fc76	add debug info	2025-05-24 15:40:07 +10:00
conrad	c0bf9d5262	group docs for find_finite() + find_nonfinite() and find_nan() + find_nonnan()	2025-05-23 12:35:23 +10:00
conradsnicta	8f74c350ce	Merge branch 'fn_find_nonnan' into 'experimental' Add `fn_find_nonnan()` See merge request conradsnicta/armadillo-code!140	2025-05-23 01:40:16 +00:00
Ryan Curtin	c9425f3158	Add `fn_find_nonnan()`	2025-05-23 01:40:16 +00:00
conrad	4e7a567967	simplifications	2025-05-22 13:29:39 +10:00
conrad	84ec3ddf9f	bump minimum gcc version to 8.1	2025-05-22 13:18:10 +10:00
conrad	0138a5b699	similar optimisation as spglue_schur	2025-05-21 10:11:43 +10:00
conrad	623eb2915d	take into account possible inf and nan	2025-05-20 22:27:38 +10:00
conrad	f845134e58	restore previous attributes due to false positives in gcc 14.2	2025-05-20 13:49:20 +10:00
conrad	014e11aa5b	remove no longer needed typedef	2025-05-20 13:20:10 +10:00
conrad	9f9a93e5e8	fix comment	2025-05-20 13:02:04 +10:00
conrad	636b8f8cf0	optimisation for empty matrices	2025-05-20 12:45:16 +10:00
conrad	d1f296eaca	simplification	2025-05-20 12:26:18 +10:00
conrad	fb5ee6fa0b	take into account possible inf and nan	2025-05-20 11:13:31 +10:00
conrad	28c3cd7e54	take into account possible inf and nan	2025-05-20 01:26:58 +10:00
conrad	69e7c73b34	take into account inf and nan	2025-05-20 01:03:12 +10:00
conrad	6fd0ff0d0f	minor layout fix	2025-05-19 12:49:13 +10:00
conrad	7e88739a41	simplification	2025-05-19 12:30:41 +10:00
conrad	0bc9d1b2d6	add explanatory comments	2025-05-19 12:22:25 +10:00
conrad	e8b7efceba	fix handling of corner-cases with nan and inf	2025-05-18 23:25:26 +10:00
conrad	7fda8d02b0	add details	2025-05-18 16:23:26 +10:00
conrad	fe3c943543	add details on unsuppported option	2025-05-15 13:30:37 +10:00
conrad	72838eb1df	better wording	2025-05-14 13:13:41 +10:00
conrad	2bf91c9159	use .is_finite()	2025-05-13 16:21:26 +10:00
conrad	c8c64e614e	conditionally undefine ARMA_IGNORE_DEPRECATED_MARKER	2025-05-13 00:29:45 +10:00
conrad	5ff113d8b1	add DOI	2025-05-12 23:30:48 +10:00
conrad	bebbf30075	deprecated syl() in favour of sylvester()	2025-05-09 10:54:26 +10:00
conrad	5b05448d69	speedup for accu(square(sp_mat.col(i)))	2025-05-09 00:41:28 +10:00
conrad	c0ef238433	deprecate stand-alone is_finite()	2025-05-07 12:45:44 +10:00
conrad	e662f86944	add TODOs	2025-05-07 12:24:39 +10:00
conrad	57da508cf3	issue compile time warning for unused outputs	2025-05-02 13:46:08 +10:00
conrad	9ac38a895c	mark as experimental	2025-04-26 21:46:56 +10:00
conrad	ed840ebb2d	use std:: prefix	2025-04-26 00:29:25 +10:00
conrad	1c41747d9c	update pdf	2025-04-25 23:56:56 +10:00
conrad	dd0aeff3f1	remove old papers	2025-04-25 23:51:16 +10:00
conrad	2186acc08f	workaround for bugs in clang 20.1	2025-04-24 15:15:37 +10:00
conrad	a2aa0f361a	rearrange	2025-04-23 16:15:18 +10:00
conrad	756cebff72	embed icon image into html	2025-04-23 14:47:45 +10:00
conrad	e767844a07	add compile time warning	2025-04-23 14:11:27 +10:00
conrad	57b7cec007	better wording	2025-04-23 13:28:54 +10:00
conrad	f3bf7bc683	add debug statement	2025-04-23 12:34:15 +10:00
conrad	a3bd59ddac	cleanup	2025-04-23 11:22:38 +10:00
conrad	82e715f2cf	clamp exponent	2025-04-23 01:27:10 +10:00
conrad	eb8ea83e11	proposed fix	2025-04-23 01:17:26 +10:00
conrad	bd912040a5	bump min version of cmake	2025-04-22 15:37:01 +10:00
conrad	f696ac0e4b	bump min version of cmake	2025-04-22 15:36:37 +10:00
conrad	616b01448f	update list of publications	2025-04-22 15:10:11 +10:00
conrad	c71edbf5e5	fixes	2025-04-22 12:51:38 +10:00
conrad	82b3972cae	add docs for ARMA_SUPERLU_64BIT_INT	2025-04-22 12:47:42 +10:00
conrad	982aaf0993	update list of supported superlu versions	2025-04-22 12:27:00 +10:00
conrad	25ba9a0ca5	state ARMA_BLAS_64BIT_INT	2025-04-22 12:24:03 +10:00
conrad	740601a431	state ARMA_BLAS_64BIT_INT	2025-04-22 12:20:41 +10:00
conrad	45423128a9	more robust handling of trace reduction	2025-04-20 22:21:35 +10:00
conrad	acd6552307	add trace reduction	2025-04-20 16:44:53 +10:00
conrad	085731e77c	increase number of terms	2025-04-20 16:31:59 +10:00
conrad	f3c1248095	allow detection of superlu 5 to 7	2025-04-19 13:33:36 +10:00
conrad	996b226da5	rename	2025-04-19 13:10:34 +10:00
conrad	3c0be74876	refactor to use superlu 7	2025-04-18 16:46:01 +10:00
conrad	aec53c9445	simplifications	2025-04-18 14:36:02 +10:00
conrad	3105fb0bd2	change error to warning	2025-04-16 00:23:18 +10:00
conrad	c2852a268d	simplifications	2025-04-11 11:51:35 +10:00
conrad	b0da624def	update docs for balance()	2025-04-09 17:22:55 +10:00
conrad	03a3efe988	update docs for balance()	2025-04-09 17:20:07 +10:00
conrad	35ba4266cd	rename variables	2025-04-09 17:15:49 +10:00
conrad	45da4eee4f	initial docs for balance()	2025-04-09 16:57:20 +10:00
conrad	90a67c55f2	move docs for affmul()	2025-04-08 16:06:36 +10:00
conrad	6346373654	prevent compiler warnings	2025-04-07 15:03:11 +10:00
conrad	7178d9715b	more debug statements	2025-04-07 14:52:47 +10:00
conrad	46fe21ae0d	improve docs for element access via better ordering	2025-04-07 14:35:49 +10:00
conrad	fc412ae23d	initial implementation of balance()	2025-04-04 23:06:36 +10:00
conrad	bc745128ef	add interfaces to xGEBAL	2025-04-04 14:46:49 +10:00
conrad	d2c7718e96	document ARMA_BLAS_64BIT_INT	2025-04-01 13:21:33 +10:00
conrad	3a0293e404	stop if both ARMA_BLAS_64BIT_INT and ARMA_USE_WRAPPER are active	2025-04-01 13:21:09 +10:00
conrad	66acb38834	expand history for new version	2025-04-01 12:55:59 +10:00
conrad	0ecd5a4c43	mark as unstable	2025-03-29 15:26:56 +10:00
conrad	0265363058	fix	2025-03-29 15:23:16 +10:00
conrad	0e091bc343	faster handling of inv(mat)*mat for tiny matrices	2025-03-29 13:03:20 +10:00