patch bump

Adds the missing argument Xst to the 'arma_ignore list'

See merge request conradsnicta/armadillo-code!142

CMakeLists.txt

@@ -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
@@ -144,7 +145,7 @@ endif()
 if(NOT APPLE)
   if(CMAKE_INSTALL_PREFIX_INITIALIZED_TO_DEFAULT)
     message(STATUS "")
-    message(STATUS "*** CMAKE_INSTALL_PREFIX was initalised by cmake to the default value of ${CMAKE_INSTALL_PREFIX}")
+    message(STATUS "*** CMAKE_INSTALL_PREFIX was initialised by cmake to the default value of ${CMAKE_INSTALL_PREFIX}")
     message(STATUS "*** CMAKE_INSTALL_PREFIX changed to /usr")
     set(CMAKE_INSTALL_PREFIX "/usr" CACHE PATH "Standard install prefix" FORCE)
   endif()
@@ -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)
@@ -618,7 +619,7 @@ install(TARGETS armadillo EXPORT ArmadilloLibraryDepends
 if(BUILD_SMOKE_TEST)
   if(HEADER_ONLY)
     ## TODO: update to link smoke_test with ARMA_LIBS when HEADER_ONLY option is enabled
-    message(STATUS "*** building smoke_test not suported when HEADER_ONLY enabled")
+    message(STATUS "*** building smoke_test not supported when HEADER_ONLY enabled")
   else()
     message(STATUS "*** configuring smoke_test")
     enable_testing()
@@ -631,7 +632,7 @@ endif()
 ## TODO: update armadillo.pc.in to link with ARMA_LIBS when HEADER_ONLY option is enabled
 
 if(HEADER_ONLY)
-  message(STATUS "*** generating pkgconfig files not suported when HEADER_ONLY enabled")
+  message(STATUS "*** generating pkgconfig files not supported when HEADER_ONLY enabled")
 else()
   message(STATUS "*** copying ${PROJECT_SOURCE_DIR}/misc/ to ${PROJECT_BINARY_DIR}/tmp/misc/")
   file(COPY ${PROJECT_SOURCE_DIR}/misc/ DESTINATION ${PROJECT_BINARY_DIR}/tmp/misc/)

NOTICE.txt

@@ -1,7 +1,7 @@
 Armadillo C++ Linear Algebra Library
-Copyright 2008-2024 Conrad Sanderson (https://conradsanderson.id.au)
+Copyright 2008-2025 Conrad Sanderson (https://conradsanderson.id.au)
 Copyright 2008-2016 National ICT Australia (NICTA)
-Copyright 2017-2024 Data61 / CSIRO
+Copyright 2017-2025 Data61 / CSIRO
 
 This product includes software developed by Conrad Sanderson (https://conradsanderson.id.au)
 This product includes software developed at National ICT Australia (NICTA)

README.md

@@ -1,9 +1,9 @@
 ### Armadillo: C++ Library for Linear Algebra & Scientific Computing  
 https://arma.sourceforge.net
 
-Copyright 2008-2024 Conrad Sanderson (https://conradsanderson.id.au)  
+Copyright 2008-2025 Conrad Sanderson (https://conradsanderson.id.au)  
 Copyright 2008-2016 National ICT Australia (NICTA)  
-Copyright 2017-2024 Data61 / CSIRO  
+Copyright 2017-2025 Data61 / CSIRO  
 
 ---
 
@@ -78,8 +78,8 @@ Please cite the following papers if you use Armadillo in your research and/or so
 Citations are useful for the continued development and maintenance of the library.
 
   * Conrad Sanderson and Ryan Curtin.  
-    Armadillo: a template-based C++ library for linear algebra.  
-    Journal of Open Source Software, Vol. 1, No. 2, pp. 26, 2016.  
+    Armadillo: An Efficient Framework for Numerical Linear Algebra.  
+    International Conference on Computer and Automation Engineering, 2025.  
   
   * Conrad Sanderson and Ryan Curtin.  
     Practical Sparse Matrices in C++ with Hybrid Storage and Template-Based Expression Optimisation.  
@@ -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`.
 
 ---

cmake_aux/Modules/ARMA_FindSuperLU.cmake

@@ -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()
+

cmake_aux/Modules/ARMA_FindSuperLU5.cmake

@@ -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()
-

examples/example1.cpp

@@ -1,7 +1,9 @@
 #include <iostream>
 #include <armadillo>
 
-using namespace std;
+using std::cout;
+using std::endl;
+
 using namespace arma;
 
 // Armadillo documentation is available at:

include/armadillo

@@ -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"

include/armadillo_bits/BaseCube_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; }
       }
     }
   

include/armadillo_bits/Base_bones.hpp

@@ -130,8 +130,8 @@ struct Base
   arma_frown("use .index_min() instead") inline elem_type min(uword& index_of_min_val) const;
   arma_frown("use .index_max() instead") inline elem_type max(uword& index_of_max_val) const;
   
-  arma_deprecated inline elem_type min(uword& row_of_min_val, uword& col_of_min_val) const;
-  arma_deprecated inline elem_type max(uword& row_of_max_val, uword& col_of_max_val) const;
+  arma_frown("use .index_min() with ind2sub() instead") inline elem_type min(uword& row_of_min_val, uword& col_of_min_val) const;
+  arma_frown("use .index_max() with ind2sub() instead") inline elem_type max(uword& row_of_max_val, uword& col_of_max_val) const;
   
   arma_warn_unused inline uword index_min() const;
   arma_warn_unused inline uword index_max() const;

include/armadillo_bits/Base_meat.hpp

@@ -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; }
         }
       }
     }

include/armadillo_bits/Col_bones.hpp

@@ -187,6 +187,8 @@ class Col<eT>::fixed : public Col<eT>
   {
   private:
   
+  using Mat<eT>::mem_local;
+  
   static constexpr bool use_extra = (fixed_n_elem > arma_config::mat_prealloc);
   
   arma_align_mem eT mem_local_extra[ (use_extra) ? fixed_n_elem : 1 ];

include/armadillo_bits/Col_meat.hpp

@@ -1345,11 +1345,11 @@ Col<eT>::fixed<fixed_n_elem>::fixed(const fill::fill_class<fill_type>&)
   {
   arma_debug_sigprint_this(this);
   
-  if(is_same_type<fill_type, fill::fill_zeros>::yes)  { (*this).zeros(); }
-  if(is_same_type<fill_type, fill::fill_ones >::yes)  { (*this).ones();  }
-  if(is_same_type<fill_type, fill::fill_eye  >::yes)  { (*this).eye();   }
-  if(is_same_type<fill_type, fill::fill_randu>::yes)  { (*this).randu(); }
-  if(is_same_type<fill_type, fill::fill_randn>::yes)  { (*this).randn(); }
+  if(is_same_type<fill_type, fill::fill_zeros>::yes)  {  (*this).zeros(); }
+  if(is_same_type<fill_type, fill::fill_ones >::yes)  {  (*this).ones();  }
+  if(is_same_type<fill_type, fill::fill_eye  >::yes)  { Mat<eT>::eye();   }
+  if(is_same_type<fill_type, fill::fill_randu>::yes)  { Mat<eT>::randu(); }
+  if(is_same_type<fill_type, fill::fill_randn>::yes)  { Mat<eT>::randn(); }
   }
 
 

include/armadillo_bits/CubeToMatOp_bones.hpp

@@ -33,8 +33,8 @@ class CubeToMatOp : public Base< typename T1::elem_type, CubeToMatOp<T1, op_type
   inline          CubeToMatOp(const T1& in_m, const uword in_aux_uword);
   inline         ~CubeToMatOp();
   
-  arma_aligned const T1&   m;            //!< the operand; must be derived from BaseCube
-  arma_aligned       uword aux_uword;    //!< auxiliary data, uword format
+  const T1&   m;            //!< the operand; must be derived from BaseCube
+        uword aux_uword;    //!< auxiliary data, uword format
   
   template<typename eT2>
   constexpr bool is_alias(const Mat<eT2>&) const { return false; }

include/armadillo_bits/Cube_bones.hpp

@@ -390,8 +390,8 @@ class Cube : public BaseCube< eT, Cube<eT> >
   arma_frown("use .index_min() instead") inline eT min(uword& index_of_min_val) const;
   arma_frown("use .index_max() instead") inline eT max(uword& index_of_max_val) const;
   
-  arma_deprecated inline eT min(uword& row_of_min_val, uword& col_of_min_val, uword& slice_of_min_val) const;
-  arma_deprecated inline eT max(uword& row_of_max_val, uword& col_of_max_val, uword& slice_of_max_val) const;
+  arma_frown("use .index_min() with ind2sub() instead") inline eT min(uword& row_of_min_val, uword& col_of_min_val, uword& slice_of_min_val) const;
+  arma_frown("use .index_max() with ind2sub() instead") inline eT max(uword& row_of_max_val, uword& col_of_max_val, uword& slice_of_max_val) const;
   
   arma_cold inline bool save(const std::string   name, const file_type type = arma_binary) const;
   arma_cold inline bool save(const hdf5_name&    spec, const file_type type = hdf5_binary) const;
@@ -485,6 +485,9 @@ class Cube<eT>::fixed : public Cube<eT>
   {
   private:
   
+  using Cube<eT>::mat_ptrs_local;
+  using Cube<eT>::mem_local;
+  
   static constexpr uword fixed_n_elem       = fixed_n_rows * fixed_n_cols * fixed_n_slices;
   static constexpr uword fixed_n_elem_slice = fixed_n_rows * fixed_n_cols;
   

include/armadillo_bits/Cube_meat.hpp

@@ -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;
   }
@@ -1019,7 +1023,7 @@ Cube<eT>::operator-=(const subview_cube<eT>& X)
 
 
 
-//! in-place element-wise cube mutiplication (using a subcube on the right-hand-side)
+//! in-place element-wise cube multiplication (using a subcube on the right-hand-side)
 template<typename eT>
 inline
 Cube<eT>&
@@ -2868,6 +2872,14 @@ Cube<eT>::Cube(const eOpCube<T1, eop_type>& X)
   
   init_cold();
   
+  if(arma_config::optimise_powexpr && is_same_type<eop_type, eop_pow>::value)
+    {
+    constexpr bool eT_non_int = is_non_integral<eT>::value;
+    
+    if(               X.aux == eT(2)   )  { eop_square::apply(*this, reinterpret_cast< const eOpCube<T1, eop_square>& >(X)); return; }
+    if(eT_non_int && (X.aux == eT(0.5)))  {   eop_sqrt::apply(*this, reinterpret_cast< const eOpCube<T1, eop_sqrt  >& >(X)); return; }
+    }
+  
   eop_type::apply(*this, X);
   }
 
@@ -2890,6 +2902,14 @@ Cube<eT>::operator=(const eOpCube<T1, eop_type>& X)
   
   init_warm(X.get_n_rows(), X.get_n_cols(), X.get_n_slices());
   
+  if(arma_config::optimise_powexpr && is_same_type<eop_type, eop_pow>::value)
+    {
+    constexpr bool eT_non_int = is_non_integral<eT>::value;
+    
+    if(               X.aux == eT(2)   )  { eop_square::apply(*this, reinterpret_cast< const eOpCube<T1, eop_square>& >(X)); return *this; }
+    if(eT_non_int && (X.aux == eT(0.5)))  {   eop_sqrt::apply(*this, reinterpret_cast< const eOpCube<T1, eop_sqrt  >& >(X)); return *this; }
+    }
+  
   eop_type::apply(*this, X);
   
   return *this;
@@ -2912,6 +2932,14 @@ Cube<eT>::operator+=(const eOpCube<T1, eop_type>& X)
   
   if(bad_alias)  { const Cube<eT> tmp(X); return (*this).operator+=(tmp); }
   
+  if(arma_config::optimise_powexpr && is_same_type<eop_type, eop_pow>::value)
+    {
+    constexpr bool eT_non_int = is_non_integral<eT>::value;
+    
+    if(               X.aux == eT(2)   )  { eop_square::apply_inplace_plus(*this, reinterpret_cast< const eOpCube<T1, eop_square>& >(X)); return *this; }
+    if(eT_non_int && (X.aux == eT(0.5)))  {   eop_sqrt::apply_inplace_plus(*this, reinterpret_cast< const eOpCube<T1, eop_sqrt  >& >(X)); return *this; }
+    }
+  
   eop_type::apply_inplace_plus(*this, X);
   
   return *this;
@@ -2934,6 +2962,14 @@ Cube<eT>::operator-=(const eOpCube<T1, eop_type>& X)
   
   if(bad_alias)  { const Cube<eT> tmp(X); return (*this).operator-=(tmp); }
   
+  if(arma_config::optimise_powexpr && is_same_type<eop_type, eop_pow>::value)
+    {
+    constexpr bool eT_non_int = is_non_integral<eT>::value;
+    
+    if(               X.aux == eT(2)   )  { eop_square::apply_inplace_minus(*this, reinterpret_cast< const eOpCube<T1, eop_square>& >(X)); return *this; }
+    if(eT_non_int && (X.aux == eT(0.5)))  {   eop_sqrt::apply_inplace_minus(*this, reinterpret_cast< const eOpCube<T1, eop_sqrt  >& >(X)); return *this; }
+    }
+  
   eop_type::apply_inplace_minus(*this, X);
   
   return *this;
@@ -2956,6 +2992,14 @@ Cube<eT>::operator%=(const eOpCube<T1, eop_type>& X)
   
   if(bad_alias)  { const Cube<eT> tmp(X); return (*this).operator%=(tmp); }
   
+  if(arma_config::optimise_powexpr && is_same_type<eop_type, eop_pow>::value)
+    {
+    constexpr bool eT_non_int = is_non_integral<eT>::value;
+    
+    if(               X.aux == eT(2)   )  { eop_square::apply_inplace_schur(*this, reinterpret_cast< const eOpCube<T1, eop_square>& >(X)); return *this; }
+    if(eT_non_int && (X.aux == eT(0.5)))  {   eop_sqrt::apply_inplace_schur(*this, reinterpret_cast< const eOpCube<T1, eop_sqrt  >& >(X)); return *this; }
+    }
+  
   eop_type::apply_inplace_schur(*this, X);
   
   return *this;
@@ -2978,6 +3022,14 @@ Cube<eT>::operator/=(const eOpCube<T1, eop_type>& X)
   
   if(bad_alias)  { const Cube<eT> tmp(X); return (*this).operator/=(tmp); }
   
+  if(arma_config::optimise_powexpr && is_same_type<eop_type, eop_pow>::value)
+    {
+    constexpr bool eT_non_int = is_non_integral<eT>::value;
+    
+    if(               X.aux == eT(2)   )  { eop_square::apply_inplace_div(*this, reinterpret_cast< const eOpCube<T1, eop_square>& >(X)); return *this; }
+    if(eT_non_int && (X.aux == eT(0.5)))  {   eop_sqrt::apply_inplace_div(*this, reinterpret_cast< const eOpCube<T1, eop_sqrt  >& >(X)); return *this; }
+    }
+  
   eop_type::apply_inplace_div(*this, X);
   
   return *this;
@@ -3448,7 +3500,7 @@ Cube<eT>::at_alt(const uword i) const
 
 
 
-//! linear element accessor (treats the cube as a vector); bounds checking not done when ARMA_NO_DEBUG is defined
+//! linear element accessor (treats the cube as a vector); bounds checking not done when ARMA_DONT_CHECK_CONFORMANCE is defined
 template<typename eT>
 arma_inline
 eT&
@@ -3461,7 +3513,7 @@ Cube<eT>::operator() (const uword i)
 
 
 
-//! linear element accessor (treats the cube as a vector); bounds checking not done when ARMA_NO_DEBUG is defined
+//! linear element accessor (treats the cube as a vector); bounds checking not done when ARMA_DONT_CHECK_CONFORMANCE is defined
 template<typename eT>
 arma_inline
 const eT&
@@ -3517,7 +3569,7 @@ Cube<eT>::at(const uword i) const
 
 
 
-//! element accessor; bounds checking not done when ARMA_NO_DEBUG is defined
+//! element accessor; bounds checking not done when ARMA_DONT_CHECK_CONFORMANCE is defined
 template<typename eT>
 arma_inline
 eT&
@@ -3537,7 +3589,7 @@ Cube<eT>::operator() (const uword in_row, const uword in_col, const uword in_sli
 
 
 
-//! element accessor; bounds checking not done when ARMA_NO_DEBUG is defined
+//! element accessor; bounds checking not done when ARMA_DONT_CHECK_CONFORMANCE is defined
 template<typename eT>
 arma_inline
 const eT&
@@ -5412,10 +5464,10 @@ Cube<eT>::fixed<fixed_n_rows, fixed_n_cols, fixed_n_slices>::fixed(const fill::f
   
   mem_setup();
   
-  if(is_same_type<fill_type, fill::fill_zeros>::yes)  { (*this).zeros(); }
-  if(is_same_type<fill_type, fill::fill_ones >::yes)  { (*this).ones();  }
-  if(is_same_type<fill_type, fill::fill_randu>::yes)  { (*this).randu(); }
-  if(is_same_type<fill_type, fill::fill_randn>::yes)  { (*this).randn(); }
+  if(is_same_type<fill_type, fill::fill_zeros>::yes)  { Cube<eT>::zeros(); }
+  if(is_same_type<fill_type, fill::fill_ones >::yes)  { Cube<eT>::ones();  }
+  if(is_same_type<fill_type, fill::fill_randu>::yes)  { Cube<eT>::randu(); }
+  if(is_same_type<fill_type, fill::fill_randn>::yes)  { Cube<eT>::randn(); }
   
   arma_static_check( (is_same_type<fill_type, fill::fill_eye>::yes), "Cube::fixed::fixed(): unsupported fill type" );
   }

include/armadillo_bits/GenCube_bones.hpp

@@ -33,9 +33,9 @@ class GenCube
   static constexpr bool use_at    = false;
   static constexpr bool is_simple = (is_same_type<gen_type, gen_ones>::value) || (is_same_type<gen_type, gen_zeros>::value);
   
-  arma_aligned const uword n_rows;
-  arma_aligned const uword n_cols;
-  arma_aligned const uword n_slices;
+  const uword n_rows;
+  const uword n_cols;
+  const uword n_slices;
   
   arma_inline  GenCube(const uword in_n_rows, const uword in_n_cols, const uword in_n_slices);
   arma_inline ~GenCube();

include/armadillo_bits/GenCube_meat.hpp

@@ -89,7 +89,7 @@ GenCube<eT, gen_type>::apply(Cube<eT>& out) const
   arma_debug_sigprint();
   
   // NOTE: we're assuming that the cube has already been set to the correct size;
-  // this is done by either the Cube contructor or operator=()
+  // this is done by either the Cube constructor or operator=()
   
        if(is_same_type<gen_type, gen_zeros>::yes) { out.zeros(); }
   else if(is_same_type<gen_type, gen_ones >::yes) { out.ones();  }

include/armadillo_bits/Gen_bones.hpp

@@ -37,8 +37,8 @@ class Gen
   static constexpr bool is_col  = T1::is_col;
   static constexpr bool is_xvec = T1::is_xvec;
   
-  arma_aligned const uword n_rows;
-  arma_aligned const uword n_cols;
+  const uword n_rows;
+  const uword n_cols;
   
   arma_inline  Gen(const uword in_n_rows, const uword in_n_cols);
   arma_inline ~Gen();

include/armadillo_bits/Gen_meat.hpp

@@ -91,7 +91,7 @@ Gen<T1, gen_type>::apply(Mat<typename T1::elem_type>& out) const
   arma_debug_sigprint();
   
   // NOTE: we're assuming that the matrix has already been set to the correct size;
-  // this is done by either the Mat contructor or operator=()
+  // this is done by either the Mat constructor or operator=()
   
        if(is_same_type<gen_type, gen_zeros>::yes) { out.zeros(); }
   else if(is_same_type<gen_type, gen_ones >::yes) { out.ones();  }

include/armadillo_bits/MapMat_meat.hpp

@@ -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 completness
-    
-    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();
         }
       }
     }

include/armadillo_bits/Mat_bones.hpp

@@ -544,8 +544,8 @@ class Mat : public Base< eT, Mat<eT> >
   arma_frown("use .index_min() instead") inline eT min(uword& index_of_min_val) const;
   arma_frown("use .index_max() instead") inline eT max(uword& index_of_max_val) const;
   
-  arma_deprecated inline eT min(uword& row_of_min_val, uword& col_of_min_val) const;
-  arma_deprecated inline eT max(uword& row_of_max_val, uword& col_of_max_val) const;
+  arma_frown("use .index_min() with ind2sub() instead") inline eT min(uword& row_of_min_val, uword& col_of_min_val) const;
+  arma_frown("use .index_max() with ind2sub() instead") inline eT max(uword& row_of_max_val, uword& col_of_max_val) const;
   
   
   arma_cold inline bool save(const std::string   name, const file_type type = arma_binary) const;
@@ -825,6 +825,8 @@ class Mat<eT>::fixed : public Mat<eT>
   {
   private:
   
+  using Mat<eT>::mem_local;
+  
   static constexpr uword fixed_n_elem = fixed_n_rows * fixed_n_cols;
   static constexpr bool  use_extra    = (fixed_n_elem > arma_config::mat_prealloc);
   

include/armadillo_bits/Mat_meat.hpp

@@ -474,9 +474,9 @@ Mat<eT>::Mat(const char* text)
   
   init( std::string(text) );
   }
-  
-  
-  
+
+
+
 //! create the matrix from a textual description
 template<typename eT>
 inline
@@ -489,8 +489,8 @@ Mat<eT>::operator=(const char* text)
   
   return *this;
   }
-  
-  
+
+
 
 //! create the matrix from a textual description
 template<typename eT>
@@ -508,9 +508,9 @@ Mat<eT>::Mat(const std::string& text)
   
   init(text);
   }
-  
-  
-  
+
+
+
 //! create the matrix from a textual description
 template<typename eT>
 inline
@@ -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;
   }
@@ -2067,7 +2071,7 @@ Mat<eT>::operator-=(const subview<eT>& X)
 
 
 
-//! in-place matrix mutiplication (using a submatrix on the right-hand-side)
+//! in-place matrix multiplication (using a submatrix on the right-hand-side)
 template<typename eT>
 inline
 Mat<eT>&
@@ -2082,7 +2086,7 @@ Mat<eT>::operator*=(const subview<eT>& X)
 
 
 
-//! in-place element-wise matrix mutiplication (using a submatrix on the right-hand-side)
+//! in-place element-wise matrix multiplication (using a submatrix on the right-hand-side)
 template<typename eT>
 inline
 Mat<eT>&
@@ -2257,7 +2261,7 @@ Mat<eT>::operator-=(const subview_cube<eT>& X)
 
 
 
-//! in-place matrix mutiplication (using a single-slice subcube on the right-hand-side)
+//! in-place matrix multiplication (using a single-slice subcube on the right-hand-side)
 template<typename eT>
 inline
 Mat<eT>&
@@ -2274,7 +2278,7 @@ Mat<eT>::operator*=(const subview_cube<eT>& X)
 
 
 
-//! in-place element-wise matrix mutiplication (using a single-slice subcube on the right-hand-side)
+//! in-place element-wise matrix multiplication (using a single-slice subcube on the right-hand-side)
 template<typename eT>
 inline
 Mat<eT>&
@@ -2383,7 +2387,7 @@ Mat<eT>::operator-=(const diagview<eT>& X)
 
 
 
-//! in-place matrix mutiplication (using a diagview on the right-hand-side)
+//! in-place matrix multiplication (using a diagview on the right-hand-side)
 template<typename eT>
 inline
 Mat<eT>&
@@ -2398,7 +2402,7 @@ Mat<eT>::operator*=(const diagview<eT>& X)
 
 
 
-//! in-place element-wise matrix mutiplication (using a diagview on the right-hand-side)
+//! in-place element-wise matrix multiplication (using a diagview on the right-hand-side)
 template<typename eT>
 inline
 Mat<eT>&
@@ -2862,6 +2866,8 @@ Mat<eT>::operator=(const SpSubview<eT>& X)
   
   if(X.n_rows == X.m.n_rows)
     {
+    arma_debug_print("access via arrays");
+    
     X.m.sync();
     
     const uword sv_col_start = X.aux_col1;
@@ -2889,6 +2895,8 @@ Mat<eT>::operator=(const SpSubview<eT>& X)
     }
   else
     {
+    arma_debug_print("access via iterators");
+    
     typename SpSubview<eT>::const_iterator it     = X.begin();
     typename SpSubview<eT>::const_iterator it_end = X.end();
     
@@ -2913,6 +2921,8 @@ Mat<eT>::operator+=(const SpSubview<eT>& X)
   
   if(X.n_rows == X.m.n_rows)
     {
+    arma_debug_print("access via arrays");
+    
     X.m.sync();
     
     const uword sv_col_start = X.aux_col1;
@@ -2940,6 +2950,8 @@ Mat<eT>::operator+=(const SpSubview<eT>& X)
     }
   else
     {
+    arma_debug_print("access via iterators");
+    
     typename SpSubview<eT>::const_iterator it     = X.begin();
     typename SpSubview<eT>::const_iterator it_end = X.end();
     
@@ -2964,6 +2976,8 @@ Mat<eT>::operator-=(const SpSubview<eT>& X)
   
   if(X.n_rows == X.m.n_rows)
     {
+    arma_debug_print("access via arrays");
+    
     X.m.sync();
     
     const uword sv_col_start = X.aux_col1;
@@ -2991,6 +3005,8 @@ Mat<eT>::operator-=(const SpSubview<eT>& X)
     }
   else
     {
+    arma_debug_print("access via iterators");
+    
     typename SpSubview<eT>::const_iterator it     = X.begin();
     typename SpSubview<eT>::const_iterator it_end = X.end();
     
@@ -5185,6 +5201,14 @@ Mat<eT>::Mat(const eOp<T1, eop_type>& X)
   
   init_cold();
   
+  if(arma_config::optimise_powexpr && is_same_type<eop_type, eop_pow>::value)
+    {
+    constexpr bool eT_non_int = is_non_integral<eT>::value;
+    
+    if(               X.aux == eT(2)   )  { eop_square::apply(*this, reinterpret_cast< const eOp<T1, eop_square>& >(X)); return; }
+    if(eT_non_int && (X.aux == eT(0.5)))  {   eop_sqrt::apply(*this, reinterpret_cast< const eOp<T1, eop_sqrt  >& >(X)); return; }
+    }
+  
   eop_type::apply(*this, X);
   }
 
@@ -5207,6 +5231,14 @@ Mat<eT>::operator=(const eOp<T1, eop_type>& X)
   
   init_warm(X.get_n_rows(), X.get_n_cols());
   
+  if(arma_config::optimise_powexpr && is_same_type<eop_type, eop_pow>::value)
+    {
+    constexpr bool eT_non_int = is_non_integral<eT>::value;
+    
+    if(               X.aux == eT(2)   )  { eop_square::apply(*this, reinterpret_cast< const eOp<T1, eop_square>& >(X)); return *this; }
+    if(eT_non_int && (X.aux == eT(0.5)))  {   eop_sqrt::apply(*this, reinterpret_cast< const eOp<T1, eop_sqrt  >& >(X)); return *this; }
+    }
+  
   eop_type::apply(*this, X);
   
   return *this;
@@ -5228,6 +5260,14 @@ Mat<eT>::operator+=(const eOp<T1, eop_type>& X)
   
   if(bad_alias)  { const Mat<eT> tmp(X); return (*this).operator+=(tmp); }
   
+  if(arma_config::optimise_powexpr && is_same_type<eop_type, eop_pow>::value)
+    {
+    constexpr bool eT_non_int = is_non_integral<eT>::value;
+    
+    if(               X.aux == eT(2)   )  { eop_square::apply_inplace_plus(*this, reinterpret_cast< const eOp<T1, eop_square>& >(X)); return *this; }
+    if(eT_non_int && (X.aux == eT(0.5)))  {   eop_sqrt::apply_inplace_plus(*this, reinterpret_cast< const eOp<T1, eop_sqrt  >& >(X)); return *this; }
+    }
+  
   eop_type::apply_inplace_plus(*this, X);
   
   return *this;
@@ -5249,6 +5289,14 @@ Mat<eT>::operator-=(const eOp<T1, eop_type>& X)
   
   if(bad_alias)  { const Mat<eT> tmp(X); return (*this).operator-=(tmp); }
   
+  if(arma_config::optimise_powexpr && is_same_type<eop_type, eop_pow>::value)
+    {
+    constexpr bool eT_non_int = is_non_integral<eT>::value;
+    
+    if(               X.aux == eT(2)   )  { eop_square::apply_inplace_minus(*this, reinterpret_cast< const eOp<T1, eop_square>& >(X)); return *this; }
+    if(eT_non_int && (X.aux == eT(0.5)))  {   eop_sqrt::apply_inplace_minus(*this, reinterpret_cast< const eOp<T1, eop_sqrt  >& >(X)); return *this; }
+    }
+  
   eop_type::apply_inplace_minus(*this, X);
   
   return *this;
@@ -5287,6 +5335,14 @@ Mat<eT>::operator%=(const eOp<T1, eop_type>& X)
   
   if(bad_alias)  { const Mat<eT> tmp(X); return (*this).operator%=(tmp); }
   
+  if(arma_config::optimise_powexpr && is_same_type<eop_type, eop_pow>::value)
+    {
+    constexpr bool eT_non_int = is_non_integral<eT>::value;
+    
+    if(               X.aux == eT(2)   )  { eop_square::apply_inplace_schur(*this, reinterpret_cast< const eOp<T1, eop_square>& >(X)); return *this; }
+    if(eT_non_int && (X.aux == eT(0.5)))  {   eop_sqrt::apply_inplace_schur(*this, reinterpret_cast< const eOp<T1, eop_sqrt  >& >(X)); return *this; }
+    }
+  
   eop_type::apply_inplace_schur(*this, X);
   
   return *this;
@@ -5308,6 +5364,14 @@ Mat<eT>::operator/=(const eOp<T1, eop_type>& X)
   
   if(bad_alias)  { const Mat<eT> tmp(X); return (*this).operator/=(tmp); }
   
+  if(arma_config::optimise_powexpr && is_same_type<eop_type, eop_pow>::value)
+    {
+    constexpr bool eT_non_int = is_non_integral<eT>::value;
+    
+    if(               X.aux == eT(2)   )  { eop_square::apply_inplace_div(*this, reinterpret_cast< const eOp<T1, eop_square>& >(X)); return *this; }
+    if(eT_non_int && (X.aux == eT(0.5)))  {   eop_sqrt::apply_inplace_div(*this, reinterpret_cast< const eOp<T1, eop_sqrt  >& >(X)); return *this; }
+    }
+  
   eop_type::apply_inplace_div(*this, X);
   
   return *this;
@@ -6391,7 +6455,7 @@ Mat<eT>::at_alt(const uword ii) const
 
 
 
-//! linear element accessor (treats the matrix as a vector); bounds checking not done when ARMA_NO_DEBUG is defined
+//! linear element accessor (treats the matrix as a vector); bounds checking not done when ARMA_DONT_CHECK_CONFORMANCE is defined
 template<typename eT>
 arma_inline
 eT&
@@ -6404,7 +6468,7 @@ Mat<eT>::operator() (const uword ii)
 
 
 
-//! linear element accessor (treats the matrix as a vector); bounds checking not done when ARMA_NO_DEBUG is defined
+//! linear element accessor (treats the matrix as a vector); bounds checking not done when ARMA_DONT_CHECK_CONFORMANCE is defined
 template<typename eT>
 arma_inline
 const eT&
@@ -6460,7 +6524,7 @@ Mat<eT>::at(const uword ii) const
 
 
 
-//! element accessor; bounds checking not done when ARMA_NO_DEBUG is defined
+//! element accessor; bounds checking not done when ARMA_DONT_CHECK_CONFORMANCE is defined
 template<typename eT>
 arma_inline
 eT&
@@ -6473,7 +6537,7 @@ Mat<eT>::operator() (const uword in_row, const uword in_col)
 
 
 
-//! element accessor; bounds checking not done when ARMA_NO_DEBUG is defined
+//! element accessor; bounds checking not done when ARMA_DONT_CHECK_CONFORMANCE is defined
 template<typename eT>
 arma_inline
 const eT&
@@ -6626,7 +6690,7 @@ Mat<eT>::is_colvec() const
 
 
 
-//! returns true if the object has the same number of non-zero rows and columnns
+//! returns true if the object has the same number of non-zero rows and columns
 template<typename eT>
 arma_inline
 bool
@@ -9565,11 +9629,11 @@ Mat<eT>::fixed<fixed_n_rows, fixed_n_cols>::fixed(const fill::fill_class<fill_ty
   {
   arma_debug_sigprint_this(this);
   
-  if(is_same_type<fill_type, fill::fill_zeros>::yes)  { (*this).zeros(); }
-  if(is_same_type<fill_type, fill::fill_ones >::yes)  { (*this).ones();  }
-  if(is_same_type<fill_type, fill::fill_eye  >::yes)  { (*this).eye();   }
-  if(is_same_type<fill_type, fill::fill_randu>::yes)  { (*this).randu(); }
-  if(is_same_type<fill_type, fill::fill_randn>::yes)  { (*this).randn(); }
+  if(is_same_type<fill_type, fill::fill_zeros>::yes)  {  (*this).zeros(); }
+  if(is_same_type<fill_type, fill::fill_ones >::yes)  {  (*this).ones();  }
+  if(is_same_type<fill_type, fill::fill_eye  >::yes)  { Mat<eT>::eye();   }
+  if(is_same_type<fill_type, fill::fill_randu>::yes)  { Mat<eT>::randu(); }
+  if(is_same_type<fill_type, fill::fill_randn>::yes)  { Mat<eT>::randn(); }
   }
 
 

include/armadillo_bits/OpCube_bones.hpp

@@ -35,11 +35,11 @@ class OpCube : public BaseCube< typename T1::elem_type, OpCube<T1, op_type> >
   inline          OpCube(const BaseCube<typename T1::elem_type, T1>& in_m, const uword in_aux_uword_a, const uword in_aux_uword_b, const uword in_aux_uword_c);
   inline         ~OpCube();
   
-  arma_aligned const T1&       m;            //!< the operand; must be derived from BaseCube
-  arma_aligned       elem_type aux;          //!< auxiliary data, using the element type as used by T1
-  arma_aligned       uword     aux_uword_a;  //!< auxiliary data, uword format
-  arma_aligned       uword     aux_uword_b;  //!< auxiliary data, uword format
-  arma_aligned       uword     aux_uword_c;  //!< auxiliary data, uword format
+  const T1&       m;            //!< the operand; must be derived from BaseCube
+        elem_type aux;          //!< auxiliary data, using the element type as used by T1
+        uword     aux_uword_a;  //!< auxiliary data, uword format
+        uword     aux_uword_b;  //!< auxiliary data, uword format
+        uword     aux_uword_c;  //!< auxiliary data, uword format
   };
 
 

include/armadillo_bits/Op_bones.hpp

@@ -61,10 +61,10 @@ class Op
   template<typename eT2>
   inline bool is_alias(const Mat<eT2>& X) const;
   
-  arma_aligned const T1&       m;            //!< the operand; must be derived from Base
-  arma_aligned       elem_type aux;          //!< auxiliary data, using the element type as used by T1
-  arma_aligned       uword     aux_uword_a;  //!< auxiliary data, uword format
-  arma_aligned       uword     aux_uword_b;  //!< auxiliary data, uword format
+  const T1&       m;            //!< the operand; must be derived from Base
+        elem_type aux;          //!< auxiliary data, using the element type as used by T1
+        uword     aux_uword_a;  //!< auxiliary data, uword format
+        uword     aux_uword_b;  //!< auxiliary data, uword format
   };
 
 

include/armadillo_bits/Proxy.hpp

@@ -168,7 +168,7 @@ struct Proxy< Mat<eT> >
   static constexpr bool is_col  = false;
   static constexpr bool is_xvec = false;
   
-  arma_aligned const Mat<eT>& Q;
+  const Mat<eT>& Q;
   
   inline explicit Proxy(const Mat<eT>& A)
     : Q(A)
@@ -215,7 +215,7 @@ struct Proxy< Col<eT> >
   static constexpr bool is_col  = true;
   static constexpr bool is_xvec = false;
   
-  arma_aligned const Col<eT>& Q;
+  const Col<eT>& Q;
   
   inline explicit Proxy(const Col<eT>& A)
     : Q(A)
@@ -262,7 +262,7 @@ struct Proxy< Row<eT> >
   static constexpr bool is_col  = false;
   static constexpr bool is_xvec = false;
   
-  arma_aligned const Row<eT>& Q;
+  const Row<eT>& Q;
   
   inline explicit Proxy(const Row<eT>& A)
     : Q(A)
@@ -309,7 +309,7 @@ struct Proxy< Gen<T1, gen_type> >
   static constexpr bool is_col  = Gen<T1, gen_type>::is_col;
   static constexpr bool is_xvec = Gen<T1, gen_type>::is_xvec;
   
-  arma_aligned const Gen<T1, gen_type>& Q;
+  const Gen<T1, gen_type>& Q;
   
   inline explicit Proxy(const Gen<T1, gen_type>& A)
     : Q(A)
@@ -356,7 +356,7 @@ struct Proxy< eOp<T1, eop_type> >
   static constexpr bool is_col  = eOp<T1, eop_type>::is_col;
   static constexpr bool is_xvec = eOp<T1, eop_type>::is_xvec;
   
-  arma_aligned const eOp<T1, eop_type>& Q;
+  const eOp<T1, eop_type>& Q;
   
   inline explicit Proxy(const eOp<T1, eop_type>& A)
     : Q(A)
@@ -403,7 +403,7 @@ struct Proxy< eGlue<T1, T2, eglue_type> >
   static constexpr bool is_col  = eGlue<T1, T2, eglue_type>::is_col;
   static constexpr bool is_xvec = eGlue<T1, T2, eglue_type>::is_xvec;
   
-  arma_aligned const eGlue<T1, T2, eglue_type>& Q;
+  const eGlue<T1, T2, eglue_type>& Q;
   
   inline explicit Proxy(const eGlue<T1, T2, eglue_type>& A)
     : Q(A)
@@ -450,7 +450,7 @@ struct Proxy< Op<T1, op_type> >
   static constexpr bool is_col  = Op<T1, op_type>::is_col;
   static constexpr bool is_xvec = Op<T1, op_type>::is_xvec;
   
-  arma_aligned const Mat<elem_type> Q;
+  const Mat<elem_type> Q;
   
   inline explicit Proxy(const Op<T1, op_type>& A)
     : Q(A)
@@ -497,7 +497,7 @@ struct Proxy< Glue<T1, T2, glue_type> >
   static constexpr bool is_col  = Glue<T1, T2, glue_type>::is_col;
   static constexpr bool is_xvec = Glue<T1, T2, glue_type>::is_xvec;
   
-  arma_aligned const Mat<elem_type> Q;
+  const Mat<elem_type> Q;
   
   inline explicit Proxy(const Glue<T1, T2, glue_type>& A)
     : Q(A)
@@ -547,9 +547,9 @@ struct Proxy< Glue<T1, T2, glue_min> >
   static constexpr bool is_col  = this_Glue_type::is_col;
   static constexpr bool is_xvec = this_Glue_type::is_xvec;
   
-  arma_aligned const this_Glue_type& Q;
-  arma_aligned const Proxy<T1>       P1;
-  arma_aligned const Proxy<T2>       P2;
+  const this_Glue_type& Q;
+  const Proxy<T1>       P1;
+  const Proxy<T2>       P2;
   
   arma_lt_comparator<elem_type> comparator;
   
@@ -605,9 +605,9 @@ struct Proxy< Glue<T1, T2, glue_max> >
   static constexpr bool is_col  = this_Glue_type::is_col;
   static constexpr bool is_xvec = this_Glue_type::is_xvec;
   
-  arma_aligned const this_Glue_type& Q;
-  arma_aligned const Proxy<T1>       P1;
-  arma_aligned const Proxy<T2>       P2;
+  const this_Glue_type& Q;
+  const Proxy<T1>       P1;
+  const Proxy<T2>       P2;
   
   arma_gt_comparator<elem_type> comparator;
   
@@ -660,7 +660,7 @@ struct Proxy< mtOp<out_eT, T1, op_type> >
   static constexpr bool is_col  = mtOp<out_eT, T1, op_type>::is_col;
   static constexpr bool is_xvec = mtOp<out_eT, T1, op_type>::is_xvec;
   
-  arma_aligned const Mat<out_eT> Q;
+  const Mat<out_eT> Q;
   
   inline explicit Proxy(const mtOp<out_eT, T1, op_type>& A)
     : Q(A)
@@ -707,7 +707,7 @@ struct Proxy< mtGlue<out_eT, T1, T2, glue_type> >
   static constexpr bool is_col  = mtGlue<out_eT, T1, T2, glue_type>::is_col;
   static constexpr bool is_xvec = mtGlue<out_eT, T1, T2, glue_type>::is_xvec;
   
-  arma_aligned const Mat<out_eT> Q;
+  const Mat<out_eT> Q;
   
   inline explicit Proxy(const mtGlue<out_eT, T1, T2, glue_type>& A)
     : Q(A)
@@ -754,7 +754,7 @@ struct Proxy< CubeToMatOp<T1, op_type> >
   static constexpr bool is_col  = CubeToMatOp<T1, op_type>::is_col;
   static constexpr bool is_xvec = CubeToMatOp<T1, op_type>::is_xvec;
   
-  arma_aligned const Mat<elem_type> Q;
+  const Mat<elem_type> Q;
   
   inline explicit Proxy(const CubeToMatOp<T1, op_type>& A)
     : Q(A)
@@ -801,8 +801,8 @@ struct Proxy< CubeToMatOp<T1, op_vectorise_cube_col> >
   static constexpr bool is_col  = true;
   static constexpr bool is_xvec = false;
   
-  arma_aligned const unwrap_cube<T1> U;
-  arma_aligned const Mat<elem_type>  Q;
+  const unwrap_cube<T1> U;
+  const Mat<elem_type>  Q;
   
   inline explicit Proxy(const CubeToMatOp<T1, op_vectorise_cube_col>& A)
     : U(A.m)
@@ -850,7 +850,7 @@ struct Proxy< SpToDOp<T1, op_type> >
   static constexpr bool is_col  = SpToDOp<T1, op_type>::is_col;
   static constexpr bool is_xvec = SpToDOp<T1, op_type>::is_xvec;
   
-  arma_aligned const Mat<elem_type> Q;
+  const Mat<elem_type> Q;
   
   inline explicit Proxy(const SpToDOp<T1, op_type>& A)
     : Q(A)
@@ -897,8 +897,8 @@ struct Proxy< SpToDOp<SpMat<eT>, op_sp_nonzeros> >
   static constexpr bool is_col  = true;
   static constexpr bool is_xvec = false;
   
-  arma_aligned const SpMat<elem_type>& R;
-  arma_aligned const Mat<elem_type>    Q;
+  const SpMat<elem_type>& R;
+  const Mat<elem_type>    Q;
   
   inline explicit Proxy(const SpToDOp<SpMat<eT>, op_sp_nonzeros>& A)
     : R(A.m)
@@ -946,7 +946,7 @@ struct Proxy< SpToDGlue<T1, T2, glue_type> >
   static constexpr bool is_col  = SpToDGlue<T1, T2, glue_type>::is_col;
   static constexpr bool is_xvec = SpToDGlue<T1, T2, glue_type>::is_xvec;
   
-  arma_aligned const Mat<elem_type> Q;
+  const Mat<elem_type> Q;
   
   inline explicit Proxy(const SpToDGlue<T1, T2, glue_type>& A)
     : Q(A)
@@ -993,7 +993,7 @@ struct Proxy< subview<eT> >
   static constexpr bool is_col  = false;
   static constexpr bool is_xvec = false;
   
-  arma_aligned const subview<eT>& Q;
+  const subview<eT>& Q;
   
   inline explicit Proxy(const subview<eT>& A)
     : Q(A)
@@ -1040,7 +1040,7 @@ struct Proxy< subview_col<eT> >
   static constexpr bool is_col  = true;
   static constexpr bool is_xvec = false;
   
-  arma_aligned const subview_col<eT>& Q;
+  const subview_col<eT>& Q;
   
   inline explicit Proxy(const subview_col<eT>& A)
     : Q(A)
@@ -1087,8 +1087,8 @@ struct Proxy< subview_cols<eT> >
   static constexpr bool is_col  = false;
   static constexpr bool is_xvec = false;
   
-  arma_aligned const subview_cols<eT>& sv;
-  arma_aligned const Mat<eT>           Q;
+  const subview_cols<eT>& sv;
+  const Mat<eT>           Q;
   
   inline explicit Proxy(const subview_cols<eT>& A)
     : sv(A)
@@ -1136,7 +1136,7 @@ struct Proxy< subview_row<eT> >
   static constexpr bool is_col  = false;
   static constexpr bool is_xvec = false;
   
-  arma_aligned const subview_row<eT>& Q;
+  const subview_row<eT>& Q;
   
   inline explicit Proxy(const subview_row<eT>& A)
     : Q(A)
@@ -1183,8 +1183,8 @@ struct Proxy< subview_elem1<eT,T1> >
   static constexpr bool is_col  = true;
   static constexpr bool is_xvec = false;
   
-  arma_aligned const subview_elem1<eT,T1>& Q;
-  arma_aligned const Proxy<T1>             R;
+  const subview_elem1<eT,T1>& Q;
+  const Proxy<T1>             R;
   
   inline explicit Proxy(const subview_elem1<eT,T1>& A)
     : Q(A)
@@ -1237,7 +1237,7 @@ struct Proxy< subview_elem2<eT,T1,T2> >
   static constexpr bool is_col  = false;
   static constexpr bool is_xvec = false;
   
-  arma_aligned const Mat<eT> Q;
+  const Mat<eT> Q;
   
   inline explicit Proxy(const subview_elem2<eT,T1,T2>& A)
     : Q(A)
@@ -1284,7 +1284,7 @@ struct Proxy< diagview<eT> >
   static constexpr bool is_col  = true;
   static constexpr bool is_xvec = false;
   
-  arma_aligned const diagview<eT>& Q;
+  const diagview<eT>& Q;
   
   inline explicit Proxy(const diagview<eT>& A)
     : Q(A)
@@ -1339,8 +1339,8 @@ struct Proxy_diagvec_mat< Op<T1, op_diagvec> >
   static constexpr bool is_col  = true;
   static constexpr bool is_xvec = false;
   
-  arma_aligned const Mat<elem_type>&     R;
-  arma_aligned const diagview<elem_type> Q;
+  const Mat<elem_type>&     R;
+  const diagview<elem_type> Q;
   
   inline explicit Proxy_diagvec_mat(const Op<T1, op_diagvec>& A)
     : R(A.m), Q( R.diag() )
@@ -1395,7 +1395,7 @@ struct Proxy_diagvec_expr< Op<T1, op_diagvec> >
   static constexpr bool is_col  = true;
   static constexpr bool is_xvec = false;
   
-  arma_aligned const Mat<elem_type> Q;
+  const Mat<elem_type> Q;
   
   inline explicit Proxy_diagvec_expr(const Op<T1, op_diagvec>& A)
     : Q(A)
@@ -1468,7 +1468,7 @@ struct Proxy< Op<T1, op_diagvec2> >
   static constexpr bool is_col  = true;
   static constexpr bool is_xvec = false;
   
-  arma_aligned const Mat<elem_type> Q;
+  const Mat<elem_type> Q;
   
   inline explicit Proxy(const Op<T1, op_diagvec2>& A)
     : Q(A)
@@ -1614,8 +1614,8 @@ struct Proxy_xtrans_vector< Op<T1, op_htrans> >
   static constexpr bool is_col  = Op<T1, op_htrans>::is_col;
   static constexpr bool is_xvec = Op<T1, op_htrans>::is_xvec;
   
-  arma_aligned const quasi_unwrap<T1> U; // avoid copy if T1 is a Row, Col or subview_col
-  arma_aligned const Mat<elem_type>   Q;
+  const quasi_unwrap<T1> U; // avoid copy if T1 is a Row, Col or subview_col
+  const Mat<elem_type>   Q;
   
   inline Proxy_xtrans_vector(const Op<T1, op_htrans>& A)
     : U(A.m)
@@ -1656,8 +1656,8 @@ struct Proxy_xtrans_vector< Op<T1, op_strans> >
   static constexpr bool is_col  = Op<T1, op_strans>::is_col;
   static constexpr bool is_xvec = Op<T1, op_strans>::is_xvec;
   
-  arma_aligned const quasi_unwrap<T1> U; // avoid copy if T1 is a Row, Col or subview_col
-  arma_aligned const Mat<elem_type>   Q;
+  const quasi_unwrap<T1> U; // avoid copy if T1 is a Row, Col or subview_col
+  const Mat<elem_type>   Q;
   
   inline Proxy_xtrans_vector(const Op<T1, op_strans>& A)
     : U(A.m)
@@ -1832,7 +1832,7 @@ struct Proxy_subview_row_htrans_cx
   static constexpr bool is_col  = true;
   static constexpr bool is_xvec = false;
   
-  arma_aligned const subview_row_htrans<eT> Q;
+  const subview_row_htrans<eT> Q;
   
   inline explicit Proxy_subview_row_htrans_cx(const Op<subview_row<eT>, op_htrans>& A)
     : Q(A.m)
@@ -1866,7 +1866,7 @@ struct Proxy_subview_row_htrans_non_cx
   static constexpr bool is_col  = true;
   static constexpr bool is_xvec = false;
   
-  arma_aligned const subview_row_strans<eT> Q;
+  const subview_row_strans<eT> Q;
   
   inline explicit Proxy_subview_row_htrans_non_cx(const Op<subview_row<eT>, op_htrans>& A)
     : Q(A.m)
@@ -1973,7 +1973,7 @@ struct Proxy< Op<subview_row<eT>, op_strans> >
   static constexpr bool is_col  = true;
   static constexpr bool is_xvec = false;
   
-  arma_aligned const subview_row_strans<eT> Q;
+  const subview_row_strans<eT> Q;
   
   inline explicit Proxy(const Op<subview_row<eT>, op_strans>& A)
     : Q(A.m)
@@ -2174,8 +2174,8 @@ struct Proxy< Op<T1, op_htrans2> >
   static constexpr bool is_col  = eOp< Op<T1, op_htrans>, eop_scalar_times>::is_col;
   static constexpr bool is_xvec = eOp< Op<T1, op_htrans>, eop_scalar_times>::is_xvec;
   
-  arma_aligned const      Op<T1, op_htrans>                     R;
-  arma_aligned const eOp< Op<T1, op_htrans>, eop_scalar_times > Q;
+  const      Op<T1, op_htrans>                     R;
+  const eOp< Op<T1, op_htrans>, eop_scalar_times > Q;
   
   inline explicit Proxy(const Op<T1, op_htrans2>& A)
     : R(A.m)
@@ -2223,7 +2223,7 @@ struct Proxy< subview_row_strans<eT> >
   static constexpr bool is_col  = true;
   static constexpr bool is_xvec = false;
   
-  arma_aligned const subview_row_strans<eT>& Q;
+  const subview_row_strans<eT>& Q;
   
   inline explicit Proxy(const subview_row_strans<eT>& A)
     : Q(A)
@@ -2270,7 +2270,7 @@ struct Proxy< subview_row_htrans<eT> >
   static constexpr bool is_col  = true;
   static constexpr bool is_xvec = false;
   
-  arma_aligned const subview_row_htrans<eT>& Q;
+  const subview_row_htrans<eT>& Q;
   
   inline explicit Proxy(const subview_row_htrans<eT>& A)
     : Q(A)
@@ -2317,7 +2317,7 @@ struct Proxy< xtrans_mat<eT, do_conj> >
   static constexpr bool is_col  = false;
   static constexpr bool is_xvec = false;
   
-  arma_aligned const Mat<eT> Q;
+  const Mat<eT> Q;
   
   inline explicit Proxy(const xtrans_mat<eT, do_conj>& A)
     : Q(A)
@@ -2364,7 +2364,7 @@ struct Proxy< xvec_htrans<eT> >
   static constexpr bool is_col  = false;
   static constexpr bool is_xvec = true;
   
-  arma_aligned const Mat<eT> Q;
+  const Mat<eT> Q;
   
   inline explicit Proxy(const xvec_htrans<eT>& A)
     : Q(A)
@@ -2419,8 +2419,8 @@ struct Proxy_vectorise_col_mat< Op<T1, op_vectorise_col> >
   static constexpr bool is_col  = true;
   static constexpr bool is_xvec = false;
   
-  arma_aligned const unwrap<T1>     U;
-  arma_aligned const Mat<elem_type> Q;
+  const unwrap<T1>     U;
+  const Mat<elem_type> Q;
   
   inline explicit Proxy_vectorise_col_mat(const Op<T1, op_vectorise_col>& A)
     : U(A.m)
@@ -2476,8 +2476,8 @@ struct Proxy_vectorise_col_expr< Op<T1, op_vectorise_col> >
   static constexpr bool is_col  = true;
   static constexpr bool is_xvec = false;
   
-  arma_aligned const Op<T1, op_vectorise_col>& Q;
-  arma_aligned const Proxy<T1>                 R;
+  const Op<T1, op_vectorise_col>& Q;
+  const Proxy<T1>                 R;
   
   inline explicit Proxy_vectorise_col_expr(const Op<T1, op_vectorise_col>& A)
     : Q(A)

include/armadillo_bits/ProxyCube.hpp

@@ -48,7 +48,7 @@ struct ProxyCube< Cube<eT> >
   static constexpr bool use_mp      = false;
   static constexpr bool has_subview = false;
   
-  arma_aligned const Cube<eT>& Q;
+  const Cube<eT>& Q;
   
   inline explicit ProxyCube(const Cube<eT>& A)
     : Q(A)
@@ -93,7 +93,7 @@ struct ProxyCube< GenCube<eT, gen_type> >
   static constexpr bool use_mp      = false;
   static constexpr bool has_subview = false;
   
-  arma_aligned const GenCube<eT, gen_type>& Q;
+  const GenCube<eT, gen_type>& Q;
   
   inline explicit ProxyCube(const GenCube<eT, gen_type>& A)
     : Q(A)
@@ -138,7 +138,7 @@ struct ProxyCube< OpCube<T1, op_type> >
   static constexpr bool use_mp      = false;
   static constexpr bool has_subview = false;
   
-  arma_aligned const Cube<elem_type> Q;
+  const Cube<elem_type> Q;
   
   inline explicit ProxyCube(const OpCube<T1, op_type>& A)
     : Q(A)
@@ -183,7 +183,7 @@ struct ProxyCube< GlueCube<T1, T2, glue_type> >
   static constexpr bool use_mp      = false;
   static constexpr bool has_subview = false;
   
-  arma_aligned const Cube<elem_type> Q;
+  const Cube<elem_type> Q;
   
   inline explicit ProxyCube(const GlueCube<T1, T2, glue_type>& A)
     : Q(A)
@@ -228,7 +228,7 @@ struct ProxyCube< subview_cube<eT> >
   static constexpr bool use_mp      = false;
   static constexpr bool has_subview = true;
   
-  arma_aligned const subview_cube<eT>& Q;
+  const subview_cube<eT>& Q;
   
   inline explicit ProxyCube(const subview_cube<eT>& A)
     : Q(A)
@@ -273,7 +273,7 @@ struct ProxyCube< subview_cube_slices<eT,T1> >
   static constexpr bool use_mp      = false;
   static constexpr bool has_subview = false;
   
-  arma_aligned const Cube<eT> Q;
+  const Cube<eT> Q;
   
   inline explicit ProxyCube(const subview_cube_slices<eT,T1>& A)
     : Q(A)
@@ -318,7 +318,7 @@ struct ProxyCube< eOpCube<T1, eop_type > >
   static constexpr bool use_mp      = eOpCube<T1, eop_type>::use_mp;
   static constexpr bool has_subview = eOpCube<T1, eop_type>::has_subview;
   
-  arma_aligned const eOpCube<T1, eop_type>& Q;
+  const eOpCube<T1, eop_type>& Q;
   
   inline explicit ProxyCube(const eOpCube<T1, eop_type>& A)
     : Q(A)
@@ -363,7 +363,7 @@ struct ProxyCube< eGlueCube<T1, T2, eglue_type > >
   static constexpr bool use_mp      = eGlueCube<T1, T2, eglue_type>::use_mp;
   static constexpr bool has_subview = eGlueCube<T1, T2, eglue_type>::has_subview;
   
-  arma_aligned const eGlueCube<T1, T2, eglue_type>& Q;
+  const eGlueCube<T1, T2, eglue_type>& Q;
   
   inline explicit ProxyCube(const eGlueCube<T1, T2, eglue_type>& A)
     : Q(A)
@@ -408,7 +408,7 @@ struct ProxyCube< mtOpCube<out_eT, T1, op_type> >
   static constexpr bool use_mp      = false;
   static constexpr bool has_subview = false;
   
-  arma_aligned const Cube<out_eT> Q;
+  const Cube<out_eT> Q;
   
   inline explicit ProxyCube(const mtOpCube<out_eT, T1, op_type>& A)
     : Q(A)
@@ -453,7 +453,7 @@ struct ProxyCube< mtGlueCube<out_eT, T1, T2, glue_type > >
   static constexpr bool use_mp      = false;
   static constexpr bool has_subview = false;
   
-  arma_aligned const Cube<out_eT> Q;
+  const Cube<out_eT> Q;
   
   inline explicit ProxyCube(const mtGlueCube<out_eT, T1, T2, glue_type>& A)
     : Q(A)

include/armadillo_bits/Row_bones.hpp

@@ -187,6 +187,8 @@ class Row<eT>::fixed : public Row<eT>
   {
   private:
   
+  using Mat<eT>::mem_local;
+  
   static constexpr bool use_extra = (fixed_n_elem > arma_config::mat_prealloc);
   
   arma_align_mem eT mem_local_extra[ (use_extra) ? fixed_n_elem : 1 ];

include/armadillo_bits/Row_meat.hpp

@@ -1345,11 +1345,11 @@ Row<eT>::fixed<fixed_n_elem>::fixed(const fill::fill_class<fill_type>&)
   {
   arma_debug_sigprint_this(this);
   
-  if(is_same_type<fill_type, fill::fill_zeros>::yes)  { (*this).zeros(); }
-  if(is_same_type<fill_type, fill::fill_ones >::yes)  { (*this).ones();  }
-  if(is_same_type<fill_type, fill::fill_eye  >::yes)  { (*this).eye();   }
-  if(is_same_type<fill_type, fill::fill_randu>::yes)  { (*this).randu(); }
-  if(is_same_type<fill_type, fill::fill_randn>::yes)  { (*this).randn(); }
+  if(is_same_type<fill_type, fill::fill_zeros>::yes)  {  (*this).zeros(); }
+  if(is_same_type<fill_type, fill::fill_ones >::yes)  {  (*this).ones();  }
+  if(is_same_type<fill_type, fill::fill_eye  >::yes)  { Mat<eT>::eye();   }
+  if(is_same_type<fill_type, fill::fill_randu>::yes)  { Mat<eT>::randu(); }
+  if(is_same_type<fill_type, fill::fill_randn>::yes)  { Mat<eT>::randn(); }
   }
 
 

include/armadillo_bits/SpBase_bones.hpp

@@ -79,8 +79,8 @@ struct SpBase
   arma_frown("use .index_min() instead") inline elem_type min(uword& index_of_min_val) const;
   arma_frown("use .index_max() instead") inline elem_type max(uword& index_of_max_val) const;
   
-  arma_deprecated inline elem_type min(uword& row_of_min_val, uword& col_of_min_val) const;
-  arma_deprecated inline elem_type max(uword& row_of_max_val, uword& col_of_max_val) const;
+  arma_frown("use .index_min() with ind2sub() instead") inline elem_type min(uword& row_of_min_val, uword& col_of_min_val) const;
+  arma_frown("use .index_max() with ind2sub() instead") inline elem_type max(uword& row_of_max_val, uword& col_of_max_val) const;
   
   arma_warn_unused inline uword index_min() const;
   arma_warn_unused inline uword index_max() const;

include/armadillo_bits/SpBase_meat.hpp

@@ -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;
       }
     }

include/armadillo_bits/SpMat_meat.hpp

@@ -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;
   }
 
@@ -3126,7 +3119,7 @@ SpMat<eT>::shed_rows(const uword in_row1, const uword in_row2)
   arma_conform_check_bounds
     (
     (in_row1 > in_row2) || (in_row2 >= n_rows),
-    "SpMat::shed_rows(): indices out of bounds or incorectly used"
+    "SpMat::shed_rows(): indices out of bounds or incorrectly used"
     );
   
   sync_csc();
@@ -3288,7 +3281,7 @@ SpMat<eT>::shed_cols(const uword in_col1, const uword in_col2)
 
 
 /**
- * Element access; acces the i'th element (works identically to the Mat accessors).
+ * Element access; access the i'th element (works identically to the Mat accessors).
  * If there is nothing at element i, 0 is returned.
  */
 
@@ -3484,7 +3477,7 @@ SpMat<eT>::is_colvec() const
 
 
 
-//! returns true if the object has the same number of non-zero rows and columnns
+//! returns true if the object has the same number of non-zero rows and columns
 template<typename eT>
 arma_inline
 bool
@@ -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;
   
@@ -5202,7 +5200,7 @@ SpMat<eT>::init(const MapMat<eT>& x)
     const uword x_index = x_entry.first;
     const eT    x_val   = x_entry.second;
     
-    // have we gone past the curent column?
+    // have we gone past the current column?
     if(x_index >= x_col_index_endp1)
       {
       x_col = x_index / x_n_rows;
@@ -5481,7 +5479,7 @@ SpMat<eT>::init_batch_add(const Mat<uword>& locs, const Mat<eT>& vals, const boo
       
       uvec sorted_indices = sort_index(abslocs); // Ascending sort.
       
-      // work out the number of unique elments 
+      // work out the number of unique elements 
       uword n_unique = 1;  // first element is unique
       
       for(uword i=1; i < sorted_indices.n_elem; ++i)
@@ -5536,7 +5534,7 @@ SpMat<eT>::init_batch_add(const Mat<uword>& locs, const Mat<eT>& vals, const boo
   
   if( (sort_locations == false) || (actually_sorted == true) )
     {
-    // work out the number of unique elments 
+    // work out the number of unique elements 
     uword n_unique = 1;  // first element is unique
     
     for(uword i=1; i < locs.n_cols; ++i)
@@ -5918,7 +5916,7 @@ SpMat<eT>::init_xform_mt(const SpBase<eT2,T1>& A, const Functor& func)
       {
       eT& t_values_i = t_values[i];
       
-      t_values_i = func(x_values[i]);   // NOTE: func() must produce a value of type eT (ie. act as a convertor between eT2 and eT)
+      t_values_i = func(x_values[i]);   // NOTE: func() must produce a value of type eT (ie. act as a converter between eT2 and eT)
       
       if(t_values_i == eT(0))  { has_zero = true; } 
       }
@@ -5936,7 +5934,7 @@ SpMat<eT>::init_xform_mt(const SpBase<eT2,T1>& A, const Functor& func)
     
     while(it != it_end)
       {
-      const eT val = func(*it);   // NOTE: func() must produce a value of type eT (ie. act as a convertor between eT2 and eT)
+      const eT val = func(*it);   // NOTE: func() must produce a value of type eT (ie. act as a converter between eT2 and eT)
       
       if(val == eT(0))  { has_zero = true; }
       

include/armadillo_bits/SpOp_bones.hpp

@@ -40,10 +40,10 @@ class SpOp : public SpBase< typename T1::elem_type, SpOp<T1, op_type> >
   
   arma_inline bool is_alias(const SpMat<elem_type>& X) const;
   
-  arma_aligned const T1&       m;            //!< the operand; must be derived from SpBase
-  arma_aligned       elem_type aux;          //!< auxiliary data, using the element type as used by T1
-  arma_aligned       uword     aux_uword_a;  //!< auxiliary data, uword format
-  arma_aligned       uword     aux_uword_b;  //!< auxiliary data, uword format
+  const T1&       m;            //!< the operand; must be derived from SpBase
+        elem_type aux;          //!< auxiliary data, using the element type as used by T1
+        uword     aux_uword_a;  //!< auxiliary data, uword format
+        uword     aux_uword_b;  //!< auxiliary data, uword format
   };
 
 

include/armadillo_bits/SpProxy.hpp

@@ -71,42 +71,42 @@ struct SpProxy< SpMat<eT> >
   typedef eT                                       elem_type;
   typedef typename get_pod_type<elem_type>::result pod_type;
   typedef SpMat<eT>                                stored_type;
-
+  
   typedef typename SpMat<eT>::const_iterator       const_iterator_type;
   typedef typename SpMat<eT>::const_row_iterator   const_row_iterator_type;
-
+  
   static constexpr bool use_iterator   = false;
   static constexpr bool Q_is_generated = false;
-
+  
   static constexpr bool is_row  = false;
   static constexpr bool is_col  = false;
   static constexpr bool is_xvec = false;
-
-  arma_aligned const SpMat<eT>& Q;
-
+  
+  const SpMat<eT>& Q;
+  
   inline explicit SpProxy(const SpMat<eT>& A)
     : Q(A)
     {
     arma_debug_sigprint();
     Q.sync();
     }
-
+  
   arma_inline uword get_n_rows()    const { return Q.n_rows;    }
   arma_inline uword get_n_cols()    const { return Q.n_cols;    }
   arma_inline uword get_n_elem()    const { return Q.n_elem;    }
   arma_inline uword get_n_nonzero() const { return Q.n_nonzero; }
-
+  
   arma_inline elem_type operator[](const uword i)                    const { return Q[i];           }
   arma_inline elem_type at        (const uword row, const uword col) const { return Q.at(row, col); }
-
+  
   arma_inline const eT*    get_values()      const { return Q.values;      }
   arma_inline const uword* get_row_indices() const { return Q.row_indices; }
   arma_inline const uword* get_col_ptrs()    const { return Q.col_ptrs;    }
-
+  
   arma_inline const_iterator_type     begin()                            const { return Q.begin();            }
   arma_inline const_iterator_type     begin_col(const uword col_num)     const { return Q.begin_col(col_num); }
   arma_inline const_row_iterator_type begin_row(const uword row_num = 0) const { return Q.begin_row(row_num); }
-
+  
   arma_inline const_iterator_type     end()                        const { return Q.end();            }
   arma_inline const_row_iterator_type end_row()                    const { return Q.end_row();        }
   arma_inline const_row_iterator_type end_row(const uword row_num) const { return Q.end_row(row_num); }
@@ -134,7 +134,7 @@ struct SpProxy< SpCol<eT> >
   static constexpr bool is_col  = true;
   static constexpr bool is_xvec = false;
   
-  arma_aligned const SpCol<eT>& Q;
+  const SpCol<eT>& Q;
   
   inline explicit SpProxy(const SpCol<eT>& A)
     : Q(A)
@@ -186,7 +186,7 @@ struct SpProxy< SpRow<eT> >
   static constexpr bool is_col  = false;
   static constexpr bool is_xvec = false;
   
-  arma_aligned const SpRow<eT>& Q;
+  const SpRow<eT>& Q;
   
   inline explicit SpProxy(const SpRow<eT>& A)
     : Q(A)
@@ -227,42 +227,42 @@ struct SpProxy< SpSubview<eT> >
   typedef eT                                           elem_type;
   typedef typename get_pod_type<elem_type>::result     pod_type;
   typedef SpSubview<eT>                                stored_type;
-
+  
   typedef typename SpSubview<eT>::const_iterator       const_iterator_type;
   typedef typename SpSubview<eT>::const_row_iterator   const_row_iterator_type;
-
+  
   static constexpr bool use_iterator   = true;
   static constexpr bool Q_is_generated = false;
-
+  
   static constexpr bool is_row  = false;
   static constexpr bool is_col  = false;
   static constexpr bool is_xvec = false;
-
-  arma_aligned const SpSubview<eT>& Q;
-
+  
+  const SpSubview<eT>& Q;
+  
   inline explicit SpProxy(const SpSubview<eT>& A)
     : Q(A)
     {
     arma_debug_sigprint();
     Q.m.sync();
     }
-
+  
   arma_inline uword get_n_rows()    const { return Q.n_rows;    }
   arma_inline uword get_n_cols()    const { return Q.n_cols;    }
   arma_inline uword get_n_elem()    const { return Q.n_elem;    }
   arma_inline uword get_n_nonzero() const { return Q.n_nonzero; }
-
+  
   arma_inline elem_type operator[](const uword i)                    const { return Q[i];           }
   arma_inline elem_type at        (const uword row, const uword col) const { return Q.at(row, col); }
-
+  
   arma_inline const eT*    get_values()      const { return Q.m.values;      }
   arma_inline const uword* get_row_indices() const { return Q.m.row_indices; }
   arma_inline const uword* get_col_ptrs()    const { return Q.m.col_ptrs;    }
-
+  
   arma_inline const_iterator_type     begin()                            const { return Q.begin();            }
   arma_inline const_iterator_type     begin_col(const uword col_num)     const { return Q.begin_col(col_num); }
   arma_inline const_row_iterator_type begin_row(const uword row_num = 0) const { return Q.begin_row(row_num); }
-
+  
   arma_inline const_iterator_type     end()                        const { return Q.end();            }
   arma_inline const_row_iterator_type end_row()                    const { return Q.end_row();        }
   arma_inline const_row_iterator_type end_row(const uword row_num) const { return Q.end_row(row_num); }
@@ -290,7 +290,7 @@ struct SpProxy< SpSubview_col<eT> >
   static constexpr bool is_col  = true;
   static constexpr bool is_xvec = false;
   
-  arma_aligned const SpSubview_col<eT>& Q;
+  const SpSubview_col<eT>& Q;
   
   inline explicit SpProxy(const SpSubview_col<eT>& A)
     : Q(A)
@@ -342,7 +342,7 @@ struct SpProxy< SpSubview_col_list<eT,T1> >
   static constexpr bool is_col  = false;
   static constexpr bool is_xvec = false;
   
-  arma_aligned const SpMat<eT> Q;
+  const SpMat<eT> Q;
   
   inline explicit SpProxy(const SpSubview_col_list<eT,T1>& A)
     : Q(A)
@@ -385,7 +385,7 @@ struct SpProxy< SpSubview_row<eT> >
   
   typedef typename SpSubview<eT>::const_iterator       const_iterator_type;
   typedef typename SpSubview<eT>::const_row_iterator   const_row_iterator_type;
-
+  
   static constexpr bool use_iterator   = true;
   static constexpr bool Q_is_generated = false;
   
@@ -393,7 +393,7 @@ struct SpProxy< SpSubview_row<eT> >
   static constexpr bool is_col  = false;
   static constexpr bool is_xvec = false;
   
-  arma_aligned const SpSubview_row<eT>& Q;
+  const SpSubview_row<eT>& Q;
   
   inline explicit SpProxy(const SpSubview_row<eT>& A)
     : Q(A)
@@ -445,7 +445,7 @@ struct SpProxy< spdiagview<eT> >
   static constexpr bool is_col  = true;
   static constexpr bool is_xvec = false;
   
-  arma_aligned const SpMat<eT> Q;
+  const SpMat<eT> Q;
   
   inline explicit SpProxy(const spdiagview<eT>& A)
     : Q(A)
@@ -497,7 +497,7 @@ struct SpProxy< SpOp<T1, spop_type> >
   static constexpr bool is_col  = SpOp<T1, spop_type>::is_col;
   static constexpr bool is_xvec = SpOp<T1, spop_type>::is_xvec;
   
-  arma_aligned const SpMat<eT> Q;
+  const SpMat<eT> Q;
   
   inline explicit SpProxy(const SpOp<T1, spop_type>& A)
     : Q(A)
@@ -549,7 +549,7 @@ struct SpProxy< SpGlue<T1, T2, spglue_type> >
   static constexpr bool is_col  = SpGlue<T1, T2, spglue_type>::is_col;
   static constexpr bool is_xvec = SpGlue<T1, T2, spglue_type>::is_xvec;
   
-  arma_aligned const SpMat<eT> Q;
+  const SpMat<eT> Q;
   
   inline explicit SpProxy(const SpGlue<T1, T2, spglue_type>& A)
     : Q(A)
@@ -600,7 +600,7 @@ struct SpProxy< mtSpOp<out_eT, T1, spop_type> >
   static constexpr bool is_col  = mtSpOp<out_eT, T1, spop_type>::is_col;
   static constexpr bool is_xvec = mtSpOp<out_eT, T1, spop_type>::is_xvec;
   
-  arma_aligned const SpMat<out_eT> Q;
+  const SpMat<out_eT> Q;
   
   inline explicit SpProxy(const mtSpOp<out_eT, T1, spop_type>& A)
     : Q(A)
@@ -651,7 +651,7 @@ struct SpProxy< mtSpGlue<out_eT, T1, T2, spglue_type> >
   static constexpr bool is_col  = mtSpGlue<out_eT, T1, T2, spglue_type>::is_col;
   static constexpr bool is_xvec = mtSpGlue<out_eT, T1, T2, spglue_type>::is_xvec;
   
-  arma_aligned const SpMat<out_eT> Q;
+  const SpMat<out_eT> Q;
   
   inline explicit SpProxy(const mtSpGlue<out_eT, T1, T2, spglue_type>& A)
     : Q(A)
@@ -702,7 +702,7 @@ struct SpProxy< mtSpReduceOp<out_eT, T1, op_type> >
   static constexpr bool is_col  = mtSpReduceOp<out_eT, T1, op_type>::is_col;
   static constexpr bool is_xvec = mtSpReduceOp<out_eT, T1, op_type>::is_xvec;
   
-  arma_aligned const SpMat<out_eT> Q;
+  const SpMat<out_eT> Q;
   
   inline explicit SpProxy(const mtSpReduceOp<out_eT, T1, op_type>& A)
     : Q(A)

include/armadillo_bits/SpSubview_col_list_meat.hpp

@@ -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;

include/armadillo_bits/SpSubview_iterators_meat.hpp

@@ -816,7 +816,7 @@ SpSubview<eT>::const_row_iterator::operator++()
         else if((*pos_ptr) == next_min_row + aux_row && col < next_min_col && (*pos_ptr) < aux_row + iterator_base::M->n_rows)
           {
           // The first element in this column is in a subsequent row that we
-          // already have another elemnt for, but the column index is less so
+          // already have another element for, but the column index is less so
           // this element will come first.
           next_min_col = col;
           next_actual_pos = col_offset + (pos_ptr - start_ptr);

include/armadillo_bits/SpSubview_meat.hpp

@@ -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();

include/armadillo_bits/SpToDOp_bones.hpp

@@ -42,10 +42,10 @@ class SpToDOp : public Base< typename T1::elem_type, SpToDOp<T1, op_type> >
   template<typename eT2>
   constexpr bool is_alias(const Mat<eT2>&) const { return false; }
   
-  arma_aligned const T1&       m;            //!< the operand; must be derived from SpBase
-  arma_aligned       elem_type aux;          //!< auxiliary data, using the element type as used by T1
-  arma_aligned       uword     aux_uword_a;  //!< auxiliary data, uword format
-  arma_aligned       uword     aux_uword_b;  //!< auxiliary data, uword format
+  const T1&       m;            //!< the operand; must be derived from SpBase
+        elem_type aux;          //!< auxiliary data, using the element type as used by T1
+        uword     aux_uword_a;  //!< auxiliary data, uword format
+        uword     aux_uword_b;  //!< auxiliary data, uword format
   };
 
 

include/armadillo_bits/SpValProxy_meat.hpp

@@ -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;

include/armadillo_bits/arma_cmath.hpp

@@ -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

include/armadillo_bits/arma_config.hpp

@@ -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

include/armadillo_bits/arma_forward.hpp

@@ -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;
@@ -130,6 +131,7 @@ class spop_strans;
 class spop_htrans;
 class spop_vectorise_row;
 class spop_vectorise_col;
+class spop_square;
 
 class spop_rel_lt_pre;
 class spop_rel_lt_post;
@@ -360,6 +362,7 @@ struct arma_zeros_indicator   : public arma_initmode_indicator<true > {};
 struct arma_nozeros_indicator : public arma_initmode_indicator<false> {};
 
 
+
 //! \addtogroup injector
 //! @{
 

include/armadillo_bits/arma_ostream_meat.hpp

@@ -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) ||

include/armadillo_bits/arma_rng.hpp

@@ -267,20 +267,14 @@ arma_rng::set_seed_random()
     {
     try
       {
-      union
-        {
-        seed_type     a;
-        unsigned char b[sizeof(seed_type)];
-        } tmp;
-      
-      tmp.a = seed_type(0);
+      char tmp[sizeof(seed_type)] = {};
       
       std::ifstream f("/dev/urandom", std::ifstream::binary);
       
-      if(f.good())  { f.read((char*)(&(tmp.b[0])), sizeof(seed_type)); }
+      if(f.good())  { f.read(&(tmp[0]), sizeof(seed_type)); }
+      
+      if(f.good())  { std::memcpy(&seed2, &(tmp[0]), sizeof(seed_type)); }
       
-      if(f.good())  { seed2 = tmp.a; }
-        
       have_seed = (seed2 != seed_type(0));
       }
     catch(...) {}
@@ -297,19 +291,17 @@ arma_rng::set_seed_random()
     
     seed3 = static_cast<seed_type>( since_epoch_usec & 0xFFFF );
     
-    union
-      {
-      uword*        a;
-      unsigned char b[sizeof(uword*)];
-      } tmp;
+    unsigned char* a = (unsigned char*)std::malloc(std::size_t(4096));
     
-    tmp.a = (uword*)malloc(sizeof(uword));
+    unsigned char  b[sizeof(unsigned char*)] = {};
     
-    if(tmp.a != nullptr)
+    if(a != nullptr)
       {
-      for(size_t i=0; i<sizeof(uword*); ++i)  { seed4 += seed_type(tmp.b[i]); }
+      std::memcpy(&(b[0]), &a, sizeof(unsigned char*));
       
-      free(tmp.a);
+      for(size_t i=0; i<sizeof(unsigned char*); ++i)  { seed4 += seed_type(b[i]); }
+      
+      std::free(a);
       }
     }
   

include/armadillo_bits/arma_version.hpp

@@ -22,9 +22,9 @@
 
 
 #define ARMA_VERSION_MAJOR 14
-#define ARMA_VERSION_MINOR 2
-#define ARMA_VERSION_PATCH 0
-#define ARMA_VERSION_NAME  "Smooth Caffeine"
+#define ARMA_VERSION_MINOR 6
+#define ARMA_VERSION_PATCH 2
+#define ARMA_VERSION_NAME  "Caffe Mocha"
 
 
 

include/armadillo_bits/arrayops_meat.hpp

@@ -239,7 +239,7 @@ arrayops::convert(out_eT* dest, const in_eT* src, const uword n_elem)
     {
     const out_eT* src2 = (const out_eT*)src;
     
-    if(dest != src2)  { arrayops::copy(dest, src2, n_elem); }
+    arrayops::copy(dest, src2, n_elem);
     
     return;
     }
@@ -304,6 +304,15 @@ inline
 void
 arrayops::convert_cx(out_eT* dest, const in_eT* src, const uword n_elem)
   {
+  if(is_same_type<out_eT,in_eT>::value)
+    {
+    const out_eT* src2 = (const out_eT*)src;
+    
+    arrayops::copy(dest, src2, n_elem);
+    
+    return;
+    }
+  
   uword j;
   
   for(j=1; j<n_elem; j+=2)
@@ -1041,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;

include/armadillo_bits/auxlib_bones.hpp

@@ -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
   

include/armadillo_bits/auxlib_meat.hpp

@@ -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
@@ -7668,6 +7738,8 @@ inline
 void_ptr
 ptr_cast(blas_int (*function)(const T*, const T*, const T*))
   {
+  // TODO: investigate replacement of union-based conversion
+  
   union converter
     {
     blas_int (*fn)(const T*, const T*, const T*);
@@ -7689,6 +7761,8 @@ inline
 void_ptr
 ptr_cast(blas_int (*function)(const std::complex<T>*, const std::complex<T>*))
   {
+  // TODO: investigate replacement of union-based conversion
+  
   union converter
     {
     blas_int (*fn)(const std::complex<T>*, const std::complex<T>*);

include/armadillo_bits/band_helper.hpp

@@ -242,7 +242,7 @@ compress(Mat<eT>& AB, const Mat<eT>& A, const uword KL, const uword KU, const bo
   // http://www.netlib.org/lapack/lug/node124.html  
   
   // for ?gbsv,  matrix AB size: 2*KL+KU+1 x N; band representation of A stored in rows KL+1 to 2*KL+KU+1  (note: fortran counts from 1)
-  // for ?gbsvx, matrix AB size:   KL+KU+1 x N; band representaiton of A stored in rows    1 to   KL+KU+1  (note: fortran counts from 1)
+  // for ?gbsvx, matrix AB size:   KL+KU+1 x N; band representation of A stored in rows    1 to   KL+KU+1  (note: fortran counts from 1)
   //
   // the +1 in the above formulas is to take into account the main diagonal
   

include/armadillo_bits/compiler_setup.hpp

@@ -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,10 +477,6 @@
 #undef major
 
 
-// optionally allow disabling of compile-time deprecation messages (not recommended)
-// NOTE: option 'ARMA_IGNORE_DEPRECATED_MARKER' will be removed
-// NOTE: disabling deprecation messages is counter-productive
-
 #if defined(ARMA_IGNORE_DEPRECATED_MARKER)
   #undef  arma_deprecated
   #define arma_deprecated

include/armadillo_bits/config.hpp

@@ -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,9 +387,12 @@
   #undef ARMA_CRIPPLED_LAPACK
 #endif
 
+// 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
 

include/armadillo_bits/config.hpp.cmake

@@ -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,9 +387,12 @@
   #undef ARMA_CRIPPLED_LAPACK
 #endif
 
+// 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
 

include/armadillo_bits/debug.hpp

@@ -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;

include/armadillo_bits/def_lapack.hpp

@@ -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
 
 
@@ -919,12 +929,18 @@ extern "C"
   void arma_fortran(arma_zhetri)(const char* uplo, const blas_int* n, blas_cxd* a, const blas_int* lda, const blas_int* ipiv, blas_cxd* work, blas_int* info, blas_len uplo_len) ARMA_NOEXCEPT;
   
   // rcond of symmetric matrix using pre-computed factorisation (real)
-  void arma_fortran(arma_ssycon)(const char* uplo, const blas_int* n, float*    a, const blas_int* lda, const blas_int* ipiv, const float*  anorm, float*  rcond, float*    work, blas_int* iwork, blas_int* info, blas_len uplo_len) ARMA_NOEXCEPT;
-  void arma_fortran(arma_dsycon)(const char* uplo, const blas_int* n, double*   a, const blas_int* lda, const blas_int* ipiv, const double* anorm, double* rcond, double*   work, blas_int* iwork, blas_int* info, blas_len uplo_len) ARMA_NOEXCEPT;
+  void arma_fortran(arma_ssycon)(const char* uplo, const blas_int* n, const float*    a, const blas_int* lda, const blas_int* ipiv, const float*  anorm, float*  rcond, float*    work, blas_int* iwork, blas_int* info, blas_len uplo_len) ARMA_NOEXCEPT;
+  void arma_fortran(arma_dsycon)(const char* uplo, const blas_int* n, const double*   a, const blas_int* lda, const blas_int* ipiv, const double* anorm, double* rcond, double*   work, blas_int* iwork, blas_int* info, blas_len uplo_len) ARMA_NOEXCEPT;
   
   // rcond of hermitian matrix using pre-computed factorisation (complex)
-  void arma_fortran(arma_checon)(const char* uplo, const blas_int* n, 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, 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;
+  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
   
@@ -1275,12 +1291,18 @@ extern "C"
   void arma_fortran(arma_zhetri)(const char* uplo, const blas_int* n, blas_cxd* a, const blas_int* lda, const blas_int* ipiv, blas_cxd* work, blas_int* info) ARMA_NOEXCEPT;
   
   // rcond of symmetric matrix using pre-computed factorisation (real)
-  void arma_fortran(arma_ssycon)(const char* uplo, const blas_int* n, float*    a, const blas_int* lda, const blas_int* ipiv, const float*  anorm,  float* rcond, float*    work, blas_int* iwork, blas_int* info) ARMA_NOEXCEPT;
-  void arma_fortran(arma_dsycon)(const char* uplo, const blas_int* n, double*   a, const blas_int* lda, const blas_int* ipiv, const double* anorm, double* rcond, double*   work, blas_int* iwork, blas_int* info) ARMA_NOEXCEPT;
+  void arma_fortran(arma_ssycon)(const char* uplo, const blas_int* n, const float*    a, const blas_int* lda, const blas_int* ipiv, const float*  anorm,  float* rcond, float*    work, blas_int* iwork, blas_int* info) ARMA_NOEXCEPT;
+  void arma_fortran(arma_dsycon)(const char* uplo, const blas_int* n, const double*   a, const blas_int* lda, const blas_int* ipiv, const double* anorm, double* rcond, double*   work, blas_int* iwork, blas_int* info) ARMA_NOEXCEPT;
   
   // rcond of hermitian matrix using pre-computed factorisation (complex)
-  void arma_fortran(arma_checon)(const char* uplo, const blas_int* n, 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, blas_cxd* a, const blas_int* lda, const blas_int* ipiv, const double* anorm, double* rcond, blas_cxd* work, blas_int* info) ARMA_NOEXCEPT;
+  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
 }

include/armadillo_bits/def_superlu.hpp

@@ -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*);

include/armadillo_bits/diagview_bones.hpp

@@ -29,7 +29,7 @@ class diagview : public Base< eT, diagview<eT> >
   typedef eT                                elem_type;
   typedef typename get_pod_type<eT>::result pod_type;
   
-  arma_aligned const Mat<eT>& m;
+  const Mat<eT>& m;
   
   static constexpr bool is_row  = false;
   static constexpr bool is_col  = true;

include/armadillo_bits/diagview_meat.hpp

@@ -586,7 +586,7 @@ diagview<eT>::extract(Mat<eT>& out, const diagview<eT>& in)
   arma_debug_sigprint();
   
   // NOTE: we're assuming that the matrix has already been set to the correct size and there is no aliasing;
-  // size setting and alias checking is done by either the Mat contructor or operator=()
+  // size setting and alias checking is done by either the Mat constructor or operator=()
   
   const Mat<eT>& in_m = in.m;
   
@@ -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];
     }
   }
 

include/armadillo_bits/diskio_meat.hpp

@@ -336,7 +336,7 @@ diskio::guess_file_type_internal(std::istream& f)
   
   // ssv_ascii has to be before csv_ascii;
   // if the data has semicolons, it suggests a CSV file with semicolon as the separating character;
-  // the semicolon may be used to allow the comma character to represent the decimal seperator (eg. 1,2345 vs 1.2345)
+  // the semicolon may be used to allow the comma character to represent the decimal separator (eg. 1,2345 vs 1.2345)
   
   if(has_semicolon && (has_bracket == false))  { return ssv_ascii; }
   
@@ -479,7 +479,7 @@ diskio::convert_token(eT& val, const std::string& token)
   //   {
   //   // std::from_chars() doesn't handle leading whitespace
   //   // std::from_chars() doesn't handle leading + sign
-  //   // std::from_chars() handles only the decimal point (.) as the decimal seperator
+  //   // std::from_chars() handles only the decimal point (.) as the decimal separator
   //   
   //   const char str0     = str[0];
   //   const bool start_ok = ((str0 != ' ') && (str0 != '\t') && (str0 != '+'));
@@ -3690,7 +3690,7 @@ diskio::load_arma_binary(SpMat<eT>& x, std::istream& f, std::string& err_msg)
         {
         arma_debug_print("detected inconsistent data while loading; re-reading integer parts as u32");
         
-        // inconstency could be due to a different uword size used during saving,
+        // inconsistency could be due to a different uword size used during saving,
         // so try loading the row_indices and col_ptrs under the assumption of 32 bit unsigned integers
         
         f.clear();
@@ -4483,7 +4483,7 @@ diskio::load_hdf5_binary(Cube<eT>& x, const hdf5_name& spec, std::string& err_ms
           return false;
           }
         
-        if(ndims == 1) { dims[1] = 1; dims[2] = 1; }  // Vector case; one row/colum, several slices
+        if(ndims == 1) { dims[1] = 1; dims[2] = 1; }  // Vector case; one row/column, several slices
         if(ndims == 2) {              dims[2] = 1; }  // Matrix case; one column, several rows/slices
         
         try { x.set_size(dims[2], dims[1], dims[0]); } catch(...) { err_msg = "not enough memory"; return false; }

include/armadillo_bits/eGlueCube_bones.hpp

@@ -27,13 +27,15 @@ class eGlueCube : public BaseCube< typename T1::elem_type, eGlueCube<T1, T2, egl
   
   typedef typename T1::elem_type                   elem_type;
   typedef typename get_pod_type<elem_type>::result pod_type;
+  typedef          ProxyCube<T1>                   proxy1_type;
+  typedef          ProxyCube<T2>                   proxy2_type;
   
   static constexpr bool use_at      = (ProxyCube<T1>::use_at      || ProxyCube<T2>::use_at     );
   static constexpr bool use_mp      = (ProxyCube<T1>::use_mp      || ProxyCube<T2>::use_mp     );
   static constexpr bool has_subview = (ProxyCube<T1>::has_subview || ProxyCube<T2>::has_subview);
   
-  arma_aligned const ProxyCube<T1> P1;
-  arma_aligned const ProxyCube<T2> P2;
+  const ProxyCube<T1> P1;
+  const ProxyCube<T2> P2;
   
   arma_inline ~eGlueCube();
   arma_inline  eGlueCube(const T1& in_A, const T2& in_B);

include/armadillo_bits/eGlue_bones.hpp

@@ -38,8 +38,8 @@ class eGlue : public Base< typename T1::elem_type, eGlue<T1, T2, eglue_type> >
   static constexpr bool is_row  = (Proxy<T1>::is_row  || Proxy<T2>::is_row );
   static constexpr bool is_xvec = (Proxy<T1>::is_xvec || Proxy<T2>::is_xvec);
   
-  arma_aligned const Proxy<T1> P1;
-  arma_aligned const Proxy<T2> P2;
+  const Proxy<T1> P1;
+  const Proxy<T2> P2;
   
   arma_inline ~eGlue();
   arma_inline  eGlue(const T1& in_A, const T2& in_B);

include/armadillo_bits/eOpCube_bones.hpp

@@ -28,16 +28,17 @@ class eOpCube : public BaseCube< typename T1::elem_type, eOpCube<T1, eop_type> >
   
   typedef typename T1::elem_type                   elem_type;
   typedef typename get_pod_type<elem_type>::result pod_type;
+  typedef          ProxyCube<T1>                   proxy_type;
   
   static constexpr bool use_at      = ProxyCube<T1>::use_at;
   static constexpr bool use_mp      = ProxyCube<T1>::use_mp || eop_type::use_mp;
   static constexpr bool has_subview = ProxyCube<T1>::has_subview;
   
-  arma_aligned const ProxyCube<T1> P;
-  arma_aligned       elem_type     aux;          //!< storage of auxiliary data, user defined format
-  arma_aligned       uword         aux_uword_a;  //!< storage of auxiliary data, uword format
-  arma_aligned       uword         aux_uword_b;  //!< storage of auxiliary data, uword format
-  arma_aligned       uword         aux_uword_c;  //!< storage of auxiliary data, uword format
+  const ProxyCube<T1> P;
+        elem_type     aux;          //!< storage of auxiliary data, user defined format
+        uword         aux_uword_a;  //!< storage of auxiliary data, uword format
+        uword         aux_uword_b;  //!< storage of auxiliary data, uword format
+        uword         aux_uword_c;  //!< storage of auxiliary data, uword format
   
   inline         ~eOpCube();
   inline explicit eOpCube(const BaseCube<typename T1::elem_type, T1>& in_m);

include/armadillo_bits/eOp_bones.hpp

@@ -38,11 +38,10 @@ class eOp : public Base< typename T1::elem_type, eOp<T1, eop_type> >
   static constexpr bool is_col  = Proxy<T1>::is_col;
   static constexpr bool is_xvec = Proxy<T1>::is_xvec;
   
-  arma_aligned const Proxy<T1> P;
-  
-  arma_aligned       elem_type aux;          //!< storage of auxiliary data, user defined format
-  arma_aligned       uword     aux_uword_a;  //!< storage of auxiliary data, uword format
-  arma_aligned       uword     aux_uword_b;  //!< storage of auxiliary data, uword format
+  const Proxy<T1> P;
+        elem_type aux;          //!< storage of auxiliary data, user defined format
+        uword     aux_uword_a;  //!< storage of auxiliary data, uword format
+        uword     aux_uword_b;  //!< storage of auxiliary data, uword format
   
   inline         ~eOp();
   inline explicit eOp(const T1& in_m);

include/armadillo_bits/eglue_core_meat.hpp

@@ -267,7 +267,7 @@ eglue_core<eglue_type>::apply(outT& out, const eGlue<T1, T2, eglue_type>& x)
   constexpr bool use_mp = (Proxy<T1>::use_mp || Proxy<T2>::use_mp) && (arma_config::openmp);
   
   // NOTE: we're assuming that the matrix has already been set to the correct size and there is no aliasing;
-  // size setting and alias checking is done by either the Mat contructor or operator=()
+  // size setting and alias checking is done by either the Mat constructor or operator=()
   
   
   eT* out_mem = out.memptr();
@@ -759,7 +759,7 @@ eglue_core<eglue_type>::apply(Cube<typename T1::elem_type>& out, const eGlueCube
   constexpr bool use_mp = (ProxyCube<T1>::use_mp || ProxyCube<T2>::use_mp) && (arma_config::openmp);
   
   // NOTE: we're assuming that the cube has already been set to the correct size and there is no aliasing;
-  // size setting and alias checking is done by either the Cube contructor or operator=()
+  // size setting and alias checking is done by either the Cube constructor or operator=()
   
   
   eT* out_mem = out.memptr();

include/armadillo_bits/eop_aux.hpp

@@ -158,7 +158,7 @@ class eop_aux
     {
     //arma_debug_sigprint();
     
-    // acording to IEEE Standard for Floating-Point Arithmetic (IEEE 754)
+    // according to IEEE Standard for Floating-Point Arithmetic (IEEE 754)
     // the mantissa length for double is 53 bits = std::numeric_limits<double>::digits
     // the mantissa length for float  is 24 bits = std::numeric_limits<float >::digits
     

include/armadillo_bits/eop_core_meat.hpp

@@ -250,7 +250,7 @@ eop_core<eop_type>::apply(outT& out, const eOp<T1, eop_type>& x)
   typedef typename T1::elem_type eT;
   
   // NOTE: we're assuming that the matrix has already been set to the correct size and there is no aliasing;
-  // size setting and alias checking is done by either the Mat contructor or operator=()
+  // size setting and alias checking is done by either the Mat constructor or operator=()
   
   const eT  k       = x.aux;
         eT* out_mem = out.memptr();
@@ -626,7 +626,7 @@ eop_core<eop_type>::apply(Cube<typename T1::elem_type>& out, const eOpCube<T1, e
   typedef typename T1::elem_type eT;
   
   // NOTE: we're assuming that the matrix has already been set to the correct size and there is no aliasing;
-  // size setting and alias checking is done by either the Mat contructor or operator=()
+  // size setting and alias checking is done by either the Mat constructor or operator=()
   
   const eT  k       = x.aux;
         eT* out_mem = out.memptr();

include/armadillo_bits/fft_engine_fftw3.hpp

@@ -49,14 +49,47 @@ class fft_engine_fftw3
   podarray<cx_type> X_work;  // for storing copy of input (can be overwritten by FFTW3)
   podarray<cx_type> Y_work;  // for storing output
   
+  inline
+  void
+  finish()
+    {
+    arma_debug_sigprint();
+    
+    if(fftw3_plan != nullptr)
+      {
+      arma_debug_print("fft_engine_fftw3::finish(): destroying plan");
+      fftw3::destroy_plan<cx_type>(fftw3_plan);
+      }
+    
+    // arma_debug_print("fft_engine_fftw3::finish(): cleanup");
+    // fftw3::cleanup<cx_type>();  // NOTE: this also removes any wisdom acquired by FFTW3 
+    }
+  
   inline
   ~fft_engine_fftw3()
     {
     arma_debug_sigprint();
     
-    if(fftw3_plan != nullptr)  { fftw3::destroy_plan<cx_type>(fftw3_plan); }
-    
-    // fftw3::cleanup<cx_type>();  // NOTE: this also removes any wisdom acquired by FFTW3 
+    #if defined(ARMA_USE_OPENMP)
+      {
+      #pragma omp critical (arma_fft_engine_fftw3)
+        {
+        (*this).finish();
+        }
+      }
+    #elif defined(ARMA_USE_STD_MUTEX)
+      {
+      std::mutex& plan_mutex = fft_engine_fftw3_aux::get_plan_mutex();
+      
+      const std::lock_guard<std::mutex> lock(plan_mutex);
+      
+      (*this).finish();
+      }
+    #else
+      {
+      (*this).finish();
+      }
+    #endif
     }
   
   inline

include/armadillo_bits/field_bones.hpp

@@ -21,13 +21,6 @@
 
 
 
-struct field_prealloc_n_elem
-  {
-  static constexpr uword val = 16;
-  };
-
-
-
 //! A lightweight 1D/2D/3D container for arbitrary objects
 //! (the objects must have a copy constructor)
 
@@ -46,8 +39,7 @@ class field
   
   private:
   
-  arma_aligned oT** mem;                                     //!< pointers to stored objects
-  arma_aligned oT*  mem_local[ field_prealloc_n_elem::val ]; //!< local storage, for small fields
+  arma_aligned oT** mem;  //!< pointers to stored objects
   
   
   public:
@@ -67,12 +59,24 @@ class field
   inline explicit field(const SizeMat&  s);
   inline explicit field(const SizeCube& s);
   
-  inline field& set_size(const uword n_obj_in);
+  inline field& set_size(const uword n_elem_in);
   inline field& set_size(const uword n_rows_in, const uword n_cols_in);
   inline field& set_size(const uword n_rows_in, const uword n_cols_in, const uword n_slices_in);
   inline field& set_size(const SizeMat&  s);
   inline field& set_size(const SizeCube& s);
   
+  inline field& reshape(const uword n_elem_in);
+  inline field& reshape(const uword n_rows_in, const uword n_cols_in);
+  inline field& reshape(const uword n_rows_in, const uword n_cols_in, const uword n_slices_in);
+  inline field& reshape(const SizeMat&  s);
+  inline field& reshape(const SizeCube& s);
+  
+  inline field& resize(const uword n_elem_in);
+  inline field& resize(const uword n_rows_in, const uword n_cols_in);
+  inline field& resize(const uword n_rows_in, const uword n_cols_in, const uword n_slices_in);
+  inline field& resize(const SizeMat&  s);
+  inline field& resize(const SizeCube& s);
+  
   inline            field(const std::vector<oT>& x);
   inline field& operator=(const std::vector<oT>& x);
   

include/armadillo_bits/field_meat.hpp

@@ -28,7 +28,7 @@ field<oT>::~field()
   
   delete_objects();
   
-  if(n_elem > field_prealloc_n_elem::val)  { delete [] mem; }
+  if(n_elem > 0)  { delete [] mem; }
   
   // try to expose buggy user code that accesses deleted objects
   mem = nullptr;
@@ -272,6 +272,173 @@ field<oT>::set_size(const SizeCube& s)
 
 
 
+template<typename oT>
+inline
+field<oT>&
+field<oT>::reshape(const uword n_elem_in)
+  {
+  arma_debug_sigprint();
+  
+  return (*this).reshape(n_elem_in, 1, 1);
+  }
+
+
+
+template<typename oT>
+inline
+field<oT>&
+field<oT>::reshape(const uword n_rows_in, const uword n_cols_in)
+  {
+  arma_debug_sigprint();
+  
+  return (*this).reshape(n_rows_in, n_cols_in, 1);
+  }
+
+
+
+template<typename oT>
+inline
+field<oT>&
+field<oT>::reshape(const uword n_rows_in, const uword n_cols_in, const uword n_slices_in)
+  {
+  arma_debug_sigprint(arma_str::format("n_rows_in: %u; n_cols_in: %u; n_slices_in: %u") % n_rows_in % n_cols_in % n_slices_in);
+  
+  if((n_rows == n_rows_in) && (n_cols == n_cols_in) && (n_slices == n_slices_in))
+    {
+    // do nothing
+    }
+  else
+  if((n_elem == 0) || ((n_rows == n_cols_in) && (n_cols == n_rows_in) && (n_slices == n_slices_in)))
+    {
+    init(n_rows_in, n_cols_in, n_slices_in);
+    }
+  else
+    {
+    field<oT> tmp(n_rows_in, n_cols_in, n_slices_in);
+    
+    const uword n_elem_to_copy = (std::min)((*this).n_elem, tmp.n_elem);
+    
+    for(uword i=0; i < n_elem_to_copy; ++i)  { tmp.at(i) = std::move((*this).at(i)); }
+    
+    (*this) = std::move(tmp);
+    }
+  
+  return *this;
+  }
+
+
+
+template<typename oT>
+inline
+field<oT>&
+field<oT>::reshape(const SizeMat& s)
+  {
+  arma_debug_sigprint();
+  
+  return (*this).reshape(s.n_rows, s.n_cols, 1);
+  }
+
+
+
+template<typename oT>
+inline
+field<oT>&
+field<oT>::reshape(const SizeCube& s)
+  {
+  arma_debug_sigprint();
+  
+  return (*this).reshape(s.n_rows, s.n_cols, s.n_slices);
+  }
+
+
+
+template<typename oT>
+inline
+field<oT>&
+field<oT>::resize(const uword n_elem_in)
+  {
+  arma_debug_sigprint();
+  
+  return (*this).resize(n_elem_in, 1, 1);
+  }
+
+
+
+template<typename oT>
+inline
+field<oT>&
+field<oT>::resize(const uword n_rows_in, const uword n_cols_in)
+  {
+  arma_debug_sigprint();
+  
+  return (*this).resize(n_rows_in, n_cols_in, 1);
+  }
+
+
+
+template<typename oT>
+inline
+field<oT>&
+field<oT>::resize(const uword n_rows_in, const uword n_cols_in, const uword n_slices_in)
+  {
+  arma_debug_sigprint(arma_str::format("n_rows_in: %u; n_cols_in: %u; n_slices_in: %u") % n_rows_in % n_cols_in % n_slices_in);
+  
+  if((n_rows == n_rows_in) && (n_cols == n_cols_in) && (n_slices == n_slices_in))
+    {
+    // do nothing
+    }
+  else
+  if(n_elem == 0)
+    {
+    (*this).set_size(n_rows_in, n_cols_in, n_slices_in);
+    }
+  else
+    {
+    // better-than-nothing implementation
+    
+    field<oT> tmp(n_rows_in, n_cols_in, n_slices_in);
+    
+    if(tmp.n_elem > 0)
+      {
+      const uword end_row   = (std::min)(n_rows_in,   n_rows  ) - 1;
+      const uword end_col   = (std::min)(n_cols_in,   n_cols  ) - 1;
+      const uword end_slice = (std::min)(n_slices_in, n_slices) - 1;
+      
+      tmp.subfield(0, 0, 0, end_row, end_col, end_slice) = (*this).subfield(0, 0, 0, end_row, end_col, end_slice);
+      }
+    
+    (*this) = std::move(tmp);
+    }
+  
+  return *this;
+  }
+
+
+
+template<typename oT>
+inline
+field<oT>&
+field<oT>::resize(const SizeMat& s)
+  {
+  arma_debug_sigprint();
+  
+  return (*this).resize(s.n_rows, s.n_cols, 1);
+  }
+
+
+
+template<typename oT>
+inline
+field<oT>&
+field<oT>::resize(const SizeCube& s)
+  {
+  arma_debug_sigprint();
+  
+  return (*this).resize(s.n_rows, s.n_cols, s.n_slices);
+  }
+
+
+
 template<typename oT>
 inline
 field<oT>::field(const std::vector<oT>& x)
@@ -279,6 +446,7 @@ field<oT>::field(const std::vector<oT>& x)
   , n_cols  (0)
   , n_slices(0)
   , n_elem  (0)
+  , mem     (nullptr)
   {
   arma_debug_sigprint_this(this);
   
@@ -312,6 +480,7 @@ field<oT>::field(const std::initializer_list<oT>& list)
   , n_cols  (0)
   , n_slices(0)
   , n_elem  (0)
+  , mem     (nullptr)
   {
   arma_debug_sigprint_this(this);
   
@@ -347,6 +516,7 @@ field<oT>::field(const std::initializer_list< std::initializer_list<oT> >& list)
   , n_cols  (0)
   , n_slices(0)
   , n_elem  (0)
+  , mem     (nullptr)
   {
   arma_debug_sigprint_this(this);
   
@@ -412,19 +582,10 @@ field<oT>::field(field<oT>&& X)
   , n_cols  (X.n_cols  )
   , n_slices(X.n_slices)
   , n_elem  (X.n_elem  )
+  , mem     (X.mem     )
   {
   arma_debug_sigprint(arma_str::format("this: %x; X: %x") % this % &X);
   
-  if(n_elem > field_prealloc_n_elem::val)
-    {
-    mem = X.mem;
-    }
-  else
-    {
-    arrayops::copy(&mem_local[0], &X.mem_local[0], n_elem);
-    mem = mem_local;
-    }
-  
   access::rw(X.n_rows  ) = 0;
   access::rw(X.n_cols  ) = 0;
   access::rw(X.n_slices) = 0;
@@ -450,15 +611,7 @@ field<oT>::operator=(field<oT>&& X)
   access::rw(n_slices) = X.n_slices;
   access::rw(n_elem  ) = X.n_elem;
   
-  if(n_elem > field_prealloc_n_elem::val)
-    {
-    mem = X.mem;
-    }
-  else
-    {
-    arrayops::copy(&mem_local[0], &X.mem_local[0], n_elem);
-    mem = mem_local;
-    }
+  mem = X.mem;
   
   access::rw(X.n_rows  ) = 0;
   access::rw(X.n_cols  ) = 0;
@@ -531,7 +684,7 @@ field<oT>::at(const uword i) const
 
 
 
-//! linear element accessor (treats the field as a vector); bounds checking not done when ARMA_NO_DEBUG is defined
+//! linear element accessor (treats the field as a vector); bounds checking not done when ARMA_DONT_CHECK_CONFORMANCE is defined
 template<typename oT>
 arma_inline
 oT&
@@ -544,7 +697,7 @@ field<oT>::operator() (const uword i)
 
 
 
-//! linear element accessor (treats the field as a vector); bounds checking not done when ARMA_NO_DEBUG is defined
+//! linear element accessor (treats the field as a vector); bounds checking not done when ARMA_DONT_CHECK_CONFORMANCE is defined
 template<typename oT>
 arma_inline
 const oT&
@@ -557,7 +710,7 @@ field<oT>::operator() (const uword i) const
 
 
 
-//! element accessor; bounds checking not done when ARMA_NO_DEBUG is defined
+//! element accessor; bounds checking not done when ARMA_DONT_CHECK_CONFORMANCE is defined
 template<typename oT>
 arma_inline
 oT&
@@ -570,7 +723,7 @@ field<oT>::operator() (const uword in_row, const uword in_col)
 
 
 
-//! element accessor; bounds checking not done when ARMA_NO_DEBUG is defined
+//! element accessor; bounds checking not done when ARMA_DONT_CHECK_CONFORMANCE is defined
 template<typename oT>
 arma_inline
 const oT&
@@ -583,7 +736,7 @@ field<oT>::operator() (const uword in_row, const uword in_col) const
 
 
 
-//! element accessor; bounds checking not done when ARMA_NO_DEBUG is defined
+//! element accessor; bounds checking not done when ARMA_DONT_CHECK_CONFORMANCE is defined
 template<typename oT>
 arma_inline
 oT&
@@ -596,7 +749,7 @@ field<oT>::operator() (const uword in_row, const uword in_col, const uword in_sl
 
 
 
-//! element accessor; bounds checking not done when ARMA_NO_DEBUG is defined
+//! element accessor; bounds checking not done when ARMA_DONT_CHECK_CONFORMANCE is defined
 template<typename oT>
 arma_inline
 const oT&
@@ -884,7 +1037,7 @@ field<oT>::rows(const uword in_row1, const uword in_row2)
   arma_conform_check_bounds
     (
     ( (in_row1 > in_row2) || (in_row2 >= n_rows) ),
-    "field::rows(): indicies out of bounds or incorrectly used"
+    "field::rows(): indices out of bounds or incorrectly used"
     );
   
   const uword sub_n_rows = in_row2 - in_row1 + 1;
@@ -907,7 +1060,7 @@ field<oT>::rows(const uword in_row1, const uword in_row2) const
   arma_conform_check_bounds
     (
     ( (in_row1 > in_row2) || (in_row2 >= n_rows) ),
-    "field::rows(): indicies out of bounds or incorrectly used"
+    "field::rows(): indices out of bounds or incorrectly used"
     );
   
   const uword sub_n_rows = in_row2 - in_row1 + 1;
@@ -930,7 +1083,7 @@ field<oT>::cols(const uword in_col1, const uword in_col2)
   arma_conform_check_bounds
     (
     ( (in_col1 > in_col2) || (in_col2 >= n_cols) ),
-    "field::cols(): indicies out of bounds or incorrectly used"
+    "field::cols(): indices out of bounds or incorrectly used"
     );
   
   const uword sub_n_cols = in_col2 - in_col1 + 1;
@@ -953,7 +1106,7 @@ field<oT>::cols(const uword in_col1, const uword in_col2) const
   arma_conform_check_bounds
     (
     ( (in_col1 > in_col2) || (in_col2 >= n_cols) ),
-    "field::cols(): indicies out of bounds or incorrectly used"
+    "field::cols(): indices out of bounds or incorrectly used"
     );
   
   const uword sub_n_cols = in_col2 - in_col1 + 1;
@@ -974,7 +1127,7 @@ field<oT>::slices(const uword in_slice1, const uword in_slice2)
   arma_conform_check_bounds
     (
     ( (in_slice1 > in_slice2) || (in_slice2 >= n_slices) ),
-    "field::slices(): indicies out of bounds or incorrectly used"
+    "field::slices(): indices out of bounds or incorrectly used"
     );
   
   const uword sub_n_slices = in_slice2 - in_slice1 + 1;
@@ -995,7 +1148,7 @@ field<oT>::slices(const uword in_slice1, const uword in_slice2) const
   arma_conform_check_bounds
     (
     ( (in_slice1 > in_slice2) || (in_slice2 >= n_slices) ),
-    "field::slices(): indicies out of bounds or incorrectly used"
+    "field::slices(): indices out of bounds or incorrectly used"
     );
   
   const uword sub_n_slices = in_slice2 - in_slice1 + 1;
@@ -1967,34 +2120,20 @@ field<oT>::init(const field<oT>& x)
   {
   arma_debug_sigprint();
   
-  if(this != &x)
+  if(this == &x)
     {
-    const uword x_n_rows   = x.n_rows;
-    const uword x_n_cols   = x.n_cols;
-    const uword x_n_slices = x.n_slices;
+    arma_debug_print("field::init(): copy omitted");
     
-    init(x_n_rows, x_n_cols, x_n_slices);
-    
-    field& t = *this;
-    
-    if(x_n_slices == 1)
-      {
-      for(uword ucol=0; ucol < x_n_cols; ++ucol)
-      for(uword urow=0; urow < x_n_rows; ++urow)
-        {
-        t.at(urow,ucol) = x.at(urow,ucol);
-        }
-      }
-    else
-      {
-      for(uword uslice=0; uslice < x_n_slices; ++uslice)
-      for(uword ucol=0;   ucol   < x_n_cols;   ++ucol  )
-      for(uword urow=0;   urow   < x_n_rows;   ++urow  )
-        {
-        t.at(urow,ucol,uslice) = x.at(urow,ucol,uslice);
-        }
-      }
+    return;
     }
+  
+  field& t = (*this);
+  
+  t.init(x.n_rows, x.n_cols, x.n_slices);
+  
+  const uword t_n_elem = t.n_elem;
+  
+  for(uword i=0; i < t_n_elem; ++i)  { t.at(i) = x.at(i); }
   }
 
 
@@ -2046,13 +2185,11 @@ field<oT>::init(const uword n_rows_in, const uword n_cols_in, const uword n_slic
     {
     delete_objects();
     
-    if(n_elem > field_prealloc_n_elem::val)  { delete [] mem; }
+    if(n_elem > 0)  { delete [] mem; }
     
-    if(n_elem_new <= field_prealloc_n_elem::val)
-      {
-      mem = (n_elem_new == 0) ? nullptr : mem_local;
-      }
-    else
+    mem = nullptr;
+    
+    if(n_elem_new > 0)
       {
       mem = new(std::nothrow) oT* [n_elem_new];
       
@@ -2079,11 +2216,7 @@ field<oT>::delete_objects()
   
   for(uword i=0; i<n_elem; ++i)
     {
-    if(mem[i] != nullptr)
-      {
-      delete mem[i];
-      mem[i] = nullptr;
-      }
+    if(mem[i] != nullptr)  { delete mem[i]; mem[i] = nullptr; }
     }
   }
 

include/armadillo_bits/fn_accu.hpp

@@ -258,20 +258,176 @@ accu(const T1& X)
   {
   arma_debug_sigprint();
   
-  const Proxy<T1> P(X);
-  
-  if(is_Mat<typename Proxy<T1>::stored_type>::value || is_subview_col<typename Proxy<T1>::stored_type>::value)
+  if((is_Mat<T1>::value) || (is_subview_col<T1>::value) || (is_Mat<typename Proxy<T1>::stored_type>::value))
     {
-    const quasi_unwrap<typename Proxy<T1>::stored_type> tmp(P.Q);
+    const quasi_unwrap<T1> U(X);
     
-    return arrayops::accumulate(tmp.M.memptr(), tmp.M.n_elem);
+    return arrayops::accumulate(U.M.memptr(), U.M.n_elem);
     }
   
+  const Proxy<T1> P(X);
+  
   return (Proxy<T1>::use_at) ? accu_proxy_at(P) : accu_proxy_linear(P);
   }
 
 
 
+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
+typename T1::elem_type
+accu(const eOp<T1,eop_square>& expr)
+  {
+  arma_debug_sigprint();
+  
+  typedef typename T1::elem_type eT;
+  
+  typedef eOp<T1,eop_square> expr_type;
+  
+  typedef typename expr_type::proxy_type::stored_type expr_P_stored_type;
+  
+  if((is_Mat<expr_P_stored_type>::value) || (is_subview_col<expr_P_stored_type>::value))
+    {
+    const quasi_unwrap<expr_P_stored_type> U(expr.P.Q);
+    
+    const eT* X_mem = U.M.memptr();
+    
+    return op_dot::direct_dot(U.M.n_elem, X_mem, X_mem);
+    }
+  
+  const Proxy<expr_type> P(expr);
+  
+  return (Proxy<expr_type>::use_at) ? accu_proxy_at(P) : accu_proxy_linear(P);
+  }
+
+
+
+template<typename T1>
+arma_warn_unused
+inline
+typename T1::elem_type
+accu(const eOp<T1,eop_pow>& expr)
+  {
+  arma_debug_sigprint();
+  
+  typedef typename T1::elem_type eT;
+  
+  typedef eOp<T1,eop_pow> expr_type;
+  
+  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(arma_config::optimise_powexpr && (expr.aux == eT(0.5)) && is_non_integral<eT>::value)
+    {
+    typedef eOp<T1,eop_sqrt> modified_expr_type;
+    
+    return accu( reinterpret_cast< const modified_expr_type& >(expr) );
+    }
+  
+  const Proxy<expr_type> P(expr);
+  
+  return (Proxy<expr_type>::use_at) ? accu_proxy_at(P) : accu_proxy_linear(P);
+  }
+
+
+
 //! explicit handling of multiply-and-accumulate
 template<typename T1, typename T2>
 arma_warn_unused
@@ -836,20 +992,83 @@ accu(const BaseCube<typename T1::elem_type,T1>& X)
   {
   arma_debug_sigprint();
   
-  const ProxyCube<T1> P(X.get_ref());
-  
-  if(is_Cube<typename ProxyCube<T1>::stored_type>::value)
+  if((is_Cube<T1>::value) || (is_Cube<typename ProxyCube<T1>::stored_type>::value))
     {
-    unwrap_cube<typename ProxyCube<T1>::stored_type> tmp(P.Q);
+    const unwrap_cube<T1> U(X.get_ref());
     
-    return arrayops::accumulate(tmp.M.memptr(), tmp.M.n_elem);
+    return arrayops::accumulate(U.M.memptr(), U.M.n_elem);
     }
   
+  const ProxyCube<T1> P(X.get_ref());
+  
   return (ProxyCube<T1>::use_at) ? accu_cube_proxy_at(P) : accu_cube_proxy_linear(P);
   }
 
 
 
+template<typename T1>
+arma_warn_unused
+inline
+typename T1::elem_type
+accu(const eOpCube<T1,eop_square>& expr)
+  {
+  arma_debug_sigprint();
+  
+  typedef typename T1::elem_type eT;
+  
+  typedef eOpCube<T1,eop_square> expr_type;
+  
+  typedef typename expr_type::proxy_type::stored_type expr_P_stored_type;
+  
+  if(is_Cube<expr_P_stored_type>::value)
+    {
+    const unwrap_cube<expr_P_stored_type> U(expr.P.Q);
+    
+    const eT* X_mem = U.M.memptr();
+    
+    return op_dot::direct_dot(U.M.n_elem, X_mem, X_mem);
+    }
+  
+  const ProxyCube<expr_type> P(expr);
+  
+  return (ProxyCube<expr_type>::use_at) ? accu_cube_proxy_at(P) : accu_cube_proxy_linear(P);
+  }
+
+
+
+template<typename T1>
+arma_warn_unused
+inline
+typename T1::elem_type
+accu(const eOpCube<T1,eop_pow>& expr)
+  {
+  arma_debug_sigprint();
+  
+  typedef typename T1::elem_type eT;
+  
+  typedef eOpCube<T1,eop_pow> expr_type;
+  
+  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(arma_config::optimise_powexpr && (expr.aux == eT(0.5)) && is_non_integral<eT>::value)
+    {
+    typedef eOpCube<T1,eop_sqrt> modified_expr_type;
+    
+    return accu( reinterpret_cast< const modified_expr_type& >(expr) );
+    }
+  
+  const ProxyCube<expr_type> P(expr);
+  
+  return (ProxyCube<expr_type>::use_at) ? accu_cube_proxy_at(P) : accu_cube_proxy_linear(P);
+  }
+
+
+
 //! explicit handling of multiply-and-accumulate (cube version)
 template<typename T1, typename T2>
 arma_warn_unused
@@ -861,8 +1080,8 @@ accu(const eGlueCube<T1,T2,eglue_schur>& expr)
   
   typedef eGlueCube<T1,T2,eglue_schur> expr_type;
   
-  typedef typename ProxyCube<T1>::stored_type P1_stored_type;
-  typedef typename ProxyCube<T2>::stored_type P2_stored_type;
+  typedef typename expr_type::proxy1_type::stored_type P1_stored_type;
+  typedef typename expr_type::proxy2_type::stored_type P2_stored_type;
   
   if(is_Cube<P1_stored_type>::value && is_Cube<P2_stored_type>::value)
     {
@@ -879,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;
+  }
+
+
+
 //
 
 
@@ -993,6 +1307,24 @@ accu(const SpGlue<T1,T2,spglue_schur>& expr)
   const SpProxy<T1> px(expr.A);
   const SpProxy<T2> py(expr.B);
   
+  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;
+  
+  if(is_SpMat<px_Q_type>::value && is_SpMat<py_Q_type>::value)
+    {
+    const unwrap_spmat<px_Q_type> UX(px.Q);
+    const unwrap_spmat<py_Q_type> UY(py.Q);
+    
+    const SpMat<eT>& X = UX.M;
+    const SpMat<eT>& Y = UY.M;
+    
+    if(&X == &Y)  { return op_dot::direct_dot(X.n_nonzero, X.values, X.values); }
+    }
+  
   typename SpProxy<T1>::const_iterator_type x_it     = px.begin();
   typename SpProxy<T1>::const_iterator_type x_it_end = px.end();
   
@@ -1020,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;
         }
       }
@@ -1058,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

include/armadillo_bits/fn_balance.hpp

@@ -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;
+  }
+
+
+
+//! @}

include/armadillo_bits/fn_conv_to.hpp

@@ -28,21 +28,44 @@ class conv_to
   {
   public:
   
-  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);
+  template<typename in_eT>
+  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_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_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_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_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_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_not_cx<in_eT>::result* junk = nullptr);
+  
+  template<typename in_eT, typename T1>
+  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);
   };
 
 
 
+template<typename out_eT>
+template<typename in_eT>
+arma_warn_unused
+inline
+out_eT
+conv_to<out_eT>::from(const in_eT& in, const typename arma_scalar_only<in_eT>::result* junk)
+  {
+  arma_debug_sigprint();
+  arma_ignore(junk);
+  
+  arma_type_check(( is_supported_elem_type<out_eT>::value == false ));
+  
+  // NOTE: this is meant only as a workaround for old user code;
+  // NOTE: it doesn't handle conversions from complex to real
+  
+  return out_eT(in);
+  }
+
+
+
 template<typename out_eT>
 template<typename in_eT, typename T1>
 arma_warn_unused

include/armadillo_bits/fn_diags_spdiags.hpp

@@ -38,7 +38,7 @@ diags(const Base<typename T1::elem_type, T1>& V_expr, const Base<sword,T2>& D_ex
   
   arma_conform_check( ((D.is_vec() == false) && (D.is_empty() == false)), "D must be a vector" );
   
-  arma_conform_check( (V.n_cols != D.n_elem), "number of colums in matrix V must match the length of vector D" );
+  arma_conform_check( (V.n_cols != D.n_elem), "number of columns in matrix V must match the length of vector D" );
   
   Mat<eT> out(n_rows, n_cols, fill::zeros);
   
@@ -93,7 +93,7 @@ spdiags(const Base<typename T1::elem_type, T1>& V_expr, const Base<sword,T2>& D_
   
   arma_conform_check( ((D.is_vec() == false) && (D.is_empty() == false)), "D must be a vector" );
   
-  arma_conform_check( (V.n_cols != D.n_elem), "number of colums in matrix V must match the length of vector D" );
+  arma_conform_check( (V.n_cols != D.n_elem), "number of columns in matrix V must match the length of vector D" );
   
   MapMat<eT> tmp(n_rows, n_cols);
   

include/armadillo_bits/fn_dot.hpp

@@ -256,7 +256,7 @@ dot
     if( &A == &B )
       {
       // We can do it directly!
-      return op_dot::direct_dot_arma(A.n_nonzero, A.values, A.values);
+      return op_dot::direct_dot(A.n_nonzero, A.values, A.values);
       }
     else
       {

include/armadillo_bits/fn_elem.hpp

@@ -52,12 +52,12 @@ real(const BaseCube<typename T1::pod_type, T1>& X)
 template<typename T1>
 arma_warn_unused
 arma_inline
-const T1&
-real(const SpBase<typename T1::pod_type,T1>& A)
+typename enable_if2< (is_arma_sparse_type<T1>::value && is_cx<typename T1::elem_type>::no), const T1& >::result
+real(const T1& X)
   {
   arma_debug_sigprint();
   
-  return A.get_ref();
+  return X;
   }
 
 
@@ -91,12 +91,12 @@ real(const BaseCube<std::complex<typename T1::pod_type>, T1>& X)
 template<typename T1>
 arma_warn_unused
 arma_inline
-const mtSpOp<typename T1::pod_type, T1, spop_real>
-real(const SpBase<std::complex<typename T1::pod_type>,T1>& A)
+typename enable_if2< (is_arma_sparse_type<T1>::value && is_cx<typename T1::elem_type>::yes), const mtSpOp<typename T1::pod_type, T1, spop_real> >::result
+real(const T1& X)
   {
   arma_debug_sigprint();
   
-  return mtSpOp<typename T1::pod_type, T1, spop_real>(A.get_ref());
+  return mtSpOp<typename T1::pod_type, T1, spop_real>(X);
   }
 
 
@@ -107,52 +107,7 @@ real(const SpBase<std::complex<typename T1::pod_type>,T1>& A)
 template<typename T1>
 arma_warn_unused
 inline
-const Gen< Mat<typename T1::pod_type>, gen_zeros >
-imag(const Base<typename T1::pod_type,T1>& X)
-  {
-  arma_debug_sigprint();
-  
-  const Proxy<T1> A(X.get_ref());
-  
-  return Gen< Mat<typename T1::pod_type>, gen_zeros>(A.get_n_rows(), A.get_n_cols());
-  }
-
-
-
-template<typename T1>
-arma_warn_unused
-inline
-const GenCube<typename T1::pod_type, gen_zeros>
-imag(const BaseCube<typename T1::pod_type,T1>& X)
-  {
-  arma_debug_sigprint();
-  
-  const ProxyCube<T1> A(X.get_ref());
-  
-  return GenCube<typename T1::pod_type, gen_zeros>(A.get_n_rows(), A.get_n_cols(), A.get_n_slices());
-  }
-
-
-
-template<typename T1>
-arma_warn_unused
-inline
-SpMat<typename T1::pod_type>
-imag(const SpBase<typename T1::pod_type,T1>& A)
-  {
-  arma_debug_sigprint();
-  
-  const SpProxy<T1> P(A.get_ref());
-  
-  return SpMat<typename T1::pod_type>(P.get_n_rows(), P.get_n_cols());
-  }
-
-
-
-template<typename T1>
-arma_warn_unused
-inline
-typename enable_if2< (is_arma_type<T1>::value && is_cx<typename T1::elem_type>::yes), const mtOp<typename T1::pod_type, T1, op_imag> >::result
+typename enable_if2< is_arma_type<T1>::value, const mtOp<typename T1::pod_type, T1, op_imag> >::result
 imag(const T1& X)
   {
   arma_debug_sigprint();
@@ -166,7 +121,7 @@ template<typename T1>
 arma_warn_unused
 inline
 const mtOpCube<typename T1::pod_type, T1, op_imag>
-imag(const BaseCube<std::complex<typename T1::pod_type>,T1>& X)
+imag(const BaseCube<typename T1::elem_type,T1>& X)
   {
   arma_debug_sigprint();
   
@@ -178,12 +133,12 @@ imag(const BaseCube<std::complex<typename T1::pod_type>,T1>& X)
 template<typename T1>
 arma_warn_unused
 arma_inline
-const mtSpOp<typename T1::pod_type, T1, spop_imag>
-imag(const SpBase<std::complex<typename T1::pod_type>,T1>& A)
+typename enable_if2< is_arma_sparse_type<T1>::value, const mtSpOp<typename T1::pod_type, T1, spop_imag> >::result
+imag(const T1& X)
   {
   arma_debug_sigprint();
   
-  return mtSpOp<typename T1::pod_type, T1, spop_imag>(A.get_ref());
+  return mtSpOp<typename T1::pod_type, T1, spop_imag>(X);
   }
 
 
@@ -480,27 +435,26 @@ abs(const BaseCube< std::complex<typename T1::pod_type>,T1>& X, const typename a
 template<typename T1>
 arma_warn_unused
 arma_inline
-const SpOp<T1, spop_abs>
-abs(const SpBase<typename T1::elem_type,T1>& X, const typename arma_not_cx<typename T1::elem_type>::result* junk = nullptr)
+typename enable_if2< (is_arma_sparse_type<T1>::value && is_cx<typename T1::elem_type>::no), const SpOp<T1, spop_abs> >::result
+abs(const T1& X)
   {
   arma_debug_sigprint();
-  arma_ignore(junk);
   
-  return SpOp<T1, spop_abs>(X.get_ref());
+  return SpOp<T1, spop_abs>(X);
   }
 
 
 
+
 template<typename T1>
 arma_warn_unused
 arma_inline
-const mtSpOp<typename T1::pod_type, T1, spop_cx_abs>
-abs(const SpBase< std::complex<typename T1::pod_type>, T1>& X, const typename arma_cx_only<typename T1::elem_type>::result* junk = nullptr)
+typename enable_if2< (is_arma_sparse_type<T1>::value && is_cx<typename T1::elem_type>::yes), const mtSpOp<typename T1::pod_type, T1, spop_cx_abs> >::result
+abs(const T1& X)
   {
   arma_debug_sigprint();
-  arma_ignore(junk);
   
-  return mtSpOp<typename T1::pod_type, T1, spop_cx_abs>(X.get_ref());
+  return mtSpOp<typename T1::pod_type, T1, spop_cx_abs>(X);
   }
 
 
@@ -568,13 +522,12 @@ arg(const BaseCube< std::complex<typename T1::pod_type>,T1>& X, const typename a
 template<typename T1>
 arma_warn_unused
 arma_inline
-const SpOp<T1, spop_arg>
-arg(const SpBase<typename T1::elem_type,T1>& X, const typename arma_not_cx<typename T1::elem_type>::result* junk = nullptr)
+typename enable_if2< (is_arma_sparse_type<T1>::value && is_cx<typename T1::elem_type>::no), const SpOp<T1, spop_arg> >::result
+arg(const T1& X)
   {
   arma_debug_sigprint();
-  arma_ignore(junk);
   
-  return SpOp<T1, spop_arg>(X.get_ref());
+  return SpOp<T1, spop_arg>(X);
   }
 
 
@@ -582,13 +535,51 @@ arg(const SpBase<typename T1::elem_type,T1>& X, const typename arma_not_cx<typen
 template<typename T1>
 arma_warn_unused
 arma_inline
-const mtSpOp<typename T1::pod_type, T1, spop_cx_arg>
-arg(const SpBase< std::complex<typename T1::pod_type>, T1>& X, const typename arma_cx_only<typename T1::elem_type>::result* junk = nullptr)
+typename enable_if2< (is_arma_sparse_type<T1>::value && is_cx<typename T1::elem_type>::yes), const mtSpOp<typename T1::pod_type, T1, spop_cx_arg> >::result
+arg(const T1& X)
   {
   arma_debug_sigprint();
-  arma_ignore(junk);
   
-  return mtSpOp<typename T1::pod_type, T1, spop_cx_arg>(X.get_ref());
+  return mtSpOp<typename T1::pod_type, T1, spop_cx_arg>(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);
   }
 
 
@@ -625,12 +616,12 @@ square(const BaseCube<typename T1::elem_type,T1>& A)
 template<typename T1>
 arma_warn_unused
 arma_inline
-const SpOp<T1, spop_square>
-square(const SpBase<typename T1::elem_type,T1>& A)
+typename enable_if2< is_arma_sparse_type<T1>::value, const SpOp<T1, spop_square> >::result
+square(const T1& X)
   {
   arma_debug_sigprint();
   
-  return SpOp<T1, spop_square>(A.get_ref());
+  return SpOp<T1, spop_square>(X);
   }
 
 
@@ -667,12 +658,12 @@ sqrt(const BaseCube<typename T1::elem_type,T1>& A)
 template<typename T1>
 arma_warn_unused
 arma_inline
-const SpOp<T1, spop_sqrt>
-sqrt(const SpBase<typename T1::elem_type,T1>& A)
+typename enable_if2< is_arma_sparse_type<T1>::value, const SpOp<T1, spop_sqrt> >::result
+sqrt(const T1& X)
   {
   arma_debug_sigprint();
   
-  return SpOp<T1, spop_sqrt>(A.get_ref());
+  return SpOp<T1, spop_sqrt>(X);
   }
 
 
@@ -709,12 +700,12 @@ cbrt(const BaseCube<typename T1::elem_type,T1>& A)
 template<typename T1>
 arma_warn_unused
 arma_inline
-typename enable_if2< is_cx<typename T1::elem_type>::no, const SpOp<T1, spop_cbrt> >::result
-cbrt(const SpBase<typename T1::elem_type,T1>& A)
+typename enable_if2< (is_arma_sparse_type<T1>::value && is_cx<typename T1::elem_type>::no), const SpOp<T1, spop_cbrt> >::result
+cbrt(const T1& X)
   {
   arma_debug_sigprint();
   
-  return SpOp<T1, spop_cbrt>(A.get_ref());
+  return SpOp<T1, spop_cbrt>(X);
   }
 
 
@@ -725,12 +716,12 @@ cbrt(const SpBase<typename T1::elem_type,T1>& A)
 template<typename T1>
 arma_warn_unused
 arma_inline
-const T1&
-conj(const Base<typename T1::pod_type,T1>& A)
+typename enable_if2< (is_arma_type<T1>::value && is_cx<typename T1::elem_type>::no), const T1& >::result
+conj(const T1& X)
   {
   arma_debug_sigprint();
   
-  return A.get_ref();
+  return X;
   }
 
 
@@ -751,12 +742,12 @@ conj(const BaseCube<typename T1::pod_type,T1>& A)
 template<typename T1>
 arma_warn_unused
 arma_inline
-const T1&
-conj(const SpBase<typename T1::pod_type,T1>& A)
+typename enable_if2< (is_arma_sparse_type<T1>::value && is_cx<typename T1::elem_type>::no), const T1& >::result
+conj(const T1& X)
   {
   arma_debug_sigprint();
   
-  return A.get_ref();
+  return X;
   }
 
 
@@ -764,12 +755,12 @@ conj(const SpBase<typename T1::pod_type,T1>& A)
 template<typename T1>
 arma_warn_unused
 arma_inline
-const eOp<T1, eop_conj>
-conj(const Base<std::complex<typename T1::pod_type>,T1>& A)
+typename enable_if2< (is_arma_type<T1>::value && is_cx<typename T1::elem_type>::yes), const eOp<T1, eop_conj> >::result
+conj(const T1& A)
   {
   arma_debug_sigprint();
   
-  return eOp<T1, eop_conj>(A.get_ref());
+  return eOp<T1, eop_conj>(A);
   }
 
 
@@ -790,12 +781,12 @@ conj(const BaseCube<std::complex<typename T1::pod_type>,T1>& A)
 template<typename T1>
 arma_warn_unused
 arma_inline
-const SpOp<T1, spop_conj>
-conj(const SpBase<std::complex<typename T1::pod_type>,T1>& A)
+typename enable_if2< (is_arma_sparse_type<T1>::value && is_cx<typename T1::elem_type>::yes), const SpOp<T1, spop_conj> >::result
+conj(const T1& X)
   {
   arma_debug_sigprint();
   
-  return SpOp<T1, spop_conj>(A.get_ref());
+  return SpOp<T1, spop_conj>(X);
   }
 
 
@@ -805,12 +796,12 @@ conj(const SpBase<std::complex<typename T1::pod_type>,T1>& A)
 template<typename T1>
 arma_warn_unused
 arma_inline
-const eOp<T1, eop_pow>
-pow(const Base<typename T1::elem_type,T1>& A, const typename T1::elem_type exponent)
+typename enable_if2< is_arma_type<T1>::value, const eOp<T1, eop_pow> >::result
+pow(const T1& A, const typename T1::elem_type exponent)
   {
   arma_debug_sigprint();
   
-  return eOp<T1, eop_pow>(A.get_ref(), exponent);
+  return eOp<T1, eop_pow>(A, exponent);
   }
 
 
@@ -833,14 +824,14 @@ pow(const BaseCube<typename T1::elem_type,T1>& A, const typename T1::elem_type e
 template<typename T1>
 arma_warn_unused
 arma_inline
-const eOp<T1, eop_pow>
-pow(const Base<typename T1::elem_type,T1>& A, const typename T1::elem_type::value_type exponent)
+typename enable_if2< (is_arma_type<T1>::value && is_cx<typename T1::elem_type>::yes), const eOp<T1, eop_pow> >::result
+pow(const T1& A, const typename T1::elem_type::value_type exponent)
   {
   arma_debug_sigprint();
   
   typedef typename T1::elem_type eT;
   
-  return eOp<T1, eop_pow>(A.get_ref(), eT(exponent));
+  return eOp<T1, eop_pow>(A, eT(exponent));
   }
 
 
@@ -849,11 +840,11 @@ template<typename T1>
 arma_warn_unused
 arma_inline
 const eOpCube<T1, eop_pow>
-pow(const BaseCube<typename T1::elem_type,T1>& A, const typename T1::elem_type::value_type exponent)
+pow(const BaseCube<std::complex<typename T1::pod_type>,T1>& A, const typename T1::elem_type::value_type exponent)
   {
   arma_debug_sigprint();
   
-  typedef typename T1::elem_type eT;
+  typedef std::complex<typename T1::pod_type> eT;
   
   return eOpCube<T1, eop_pow>(A.get_ref(), eT(exponent));
   }
@@ -892,12 +883,12 @@ floor(const BaseCube<typename T1::elem_type,T1>& A)
 template<typename T1>
 arma_warn_unused
 arma_inline
-const SpOp<T1, spop_floor>
-floor(const SpBase<typename T1::elem_type,T1>& X)
+typename enable_if2< is_arma_sparse_type<T1>::value, const SpOp<T1, spop_floor> >::result
+floor(const T1& X)
   {
   arma_debug_sigprint();
   
-  return SpOp<T1, spop_floor>(X.get_ref());
+  return SpOp<T1, spop_floor>(X);
   }
 
 
@@ -934,12 +925,12 @@ ceil(const BaseCube<typename T1::elem_type,T1>& A)
 template<typename T1>
 arma_warn_unused
 arma_inline
-const SpOp<T1, spop_ceil>
-ceil(const SpBase<typename T1::elem_type,T1>& X)
+typename enable_if2< is_arma_sparse_type<T1>::value, const SpOp<T1, spop_ceil> >::result
+ceil(const T1& X)
   {
   arma_debug_sigprint();
   
-  return SpOp<T1, spop_ceil>(X.get_ref());
+  return SpOp<T1, spop_ceil>(X);
   }
 
 
@@ -976,12 +967,12 @@ round(const BaseCube<typename T1::elem_type,T1>& A)
 template<typename T1>
 arma_warn_unused
 arma_inline
-const SpOp<T1, spop_round>
-round(const SpBase<typename T1::elem_type,T1>& X)
+typename enable_if2< is_arma_sparse_type<T1>::value, const SpOp<T1, spop_round> >::result
+round(const T1& X)
   {
   arma_debug_sigprint();
   
-  return SpOp<T1, spop_round>(X.get_ref());
+  return SpOp<T1, spop_round>(X);
   }
 
 
@@ -1018,12 +1009,12 @@ trunc(const BaseCube<typename T1::elem_type,T1>& A)
 template<typename T1>
 arma_warn_unused
 arma_inline
-const SpOp<T1, spop_trunc>
-trunc(const SpBase<typename T1::elem_type,T1>& X)
+typename enable_if2< is_arma_sparse_type<T1>::value, const SpOp<T1, spop_trunc> >::result
+trunc(const T1& X)
   {
   arma_debug_sigprint();
   
-  return SpOp<T1, spop_trunc>(X.get_ref());
+  return SpOp<T1, spop_trunc>(X);
   }
 
 
@@ -1073,12 +1064,12 @@ sign(const BaseCube<typename T1::elem_type,T1>& A)
 template<typename T1>
 arma_warn_unused
 arma_inline
-const SpOp<T1, spop_sign>
-sign(const SpBase<typename T1::elem_type,T1>& X)
+typename enable_if2< is_arma_sparse_type<T1>::value, const SpOp<T1, spop_sign> >::result
+sign(const T1& X)
   {
   arma_debug_sigprint();
   
-  return SpOp<T1, spop_sign>(X.get_ref());
+  return SpOp<T1, spop_sign>(X);
   }
 
 

include/armadillo_bits/fn_find.hpp

@@ -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;
+  }
+
+
+
 //! @}

include/armadillo_bits/fn_interp1.hpp

@@ -321,7 +321,7 @@ interp1
     else if(c1 == 'l')  { sig = 20; }  // linear
     else
       {
-      if( (c1 == '*') && (c2 == 'n') )  { sig = 11; }  // nearest neighour, assume monotonic increase in X and XI
+      if( (c1 == '*') && (c2 == 'n') )  { sig = 11; }  // nearest neighbour, assume monotonic increase in X and XI
       if( (c1 == '*') && (c2 == 'l') )  { sig = 21; }  // linear, assume monotonic increase in X and XI
       }
     }

include/armadillo_bits/fn_misc.hpp

@@ -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;
       }

include/armadillo_bits/fn_omit.hpp

@@ -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);
+  }
+
+
+
+//! @}

include/armadillo_bits/fn_powext.hpp

@@ -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

include/armadillo_bits/fn_randperm.hpp

@@ -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;
     }
   

include/armadillo_bits/fn_reshape.hpp

@@ -135,4 +135,53 @@ reshape(const SpBase<typename T1::elem_type, T1>& X, const SizeMat& s)
 
 
 
+//
+
+
+
+template<typename oT>
+arma_warn_unused
+inline
+field<oT>
+reshape(const field<oT>& A, const uword new_n_rows, const uword new_n_cols, const uword new_n_slices = uword(1))
+  {
+  arma_debug_sigprint();
+  
+  field<oT> B(new_n_rows, new_n_cols, new_n_slices);
+  
+  const uword n_elem_to_copy = (std::min)(A.n_elem, B.n_elem);
+  
+  for(uword i=0; i < n_elem_to_copy; ++i)  { B.at(i) = A.at(i); }
+  
+  return B;
+  }
+
+
+
+template<typename oT>
+arma_warn_unused
+inline
+field<oT>
+reshape(const field<oT>& A, const SizeMat& s)
+  {
+  arma_debug_sigprint();
+  
+  return reshape(A, s.n_rows, s.n_cols);
+  }
+
+
+
+template<typename oT>
+arma_warn_unused
+inline
+field<oT>
+reshape(const field<oT>& A, const SizeCube& s)
+  {
+  arma_debug_sigprint();
+  
+  return reshape(A, s.n_rows, s.n_cols, s.n_slices);
+  }
+
+
+
 //! @}

include/armadillo_bits/fn_resize.hpp

@@ -99,4 +99,56 @@ resize(const SpBase<typename T1::elem_type, T1>& X, const SizeMat& s)
 
 
 
+template<typename oT>
+arma_warn_unused
+inline
+field<oT>
+resize(const field<oT>& A, const uword in_n_rows, const uword in_n_cols, const uword in_n_slices = uword(1))
+  {
+  arma_debug_sigprint();
+  
+  // better-than-nothing implementation
+  
+  field<oT> B(in_n_rows, in_n_cols, in_n_slices);
+  
+  if((B.n_elem > 0) && (A.n_elem > 0))
+    {
+    const uword end_row   = (std::min)(in_n_rows,   A.n_rows  ) - 1;
+    const uword end_col   = (std::min)(in_n_cols,   A.n_cols  ) - 1;
+    const uword end_slice = (std::min)(in_n_slices, A.n_slices) - 1;
+    
+    B.subfield(0, 0, 0, end_row, end_col, end_slice) = A.subfield(0, 0, 0, end_row, end_col, end_slice);
+    }
+  
+  return B;
+  }
+
+
+
+template<typename oT>
+arma_warn_unused
+inline
+field<oT>
+resize(const field<oT>& A, const SizeMat& s)
+  {
+  arma_debug_sigprint();
+  
+  return resize(A, s.n_rows, s.n_cols);
+  }
+
+
+
+template<typename oT>
+arma_warn_unused
+inline
+field<oT>
+resize(const field<oT>& A, const SizeCube& s)
+  {
+  arma_debug_sigprint();
+  
+  return resize(A, s.n_rows, s.n_cols, s.n_slices);
+  }
+
+
+
 //! @}

include/armadillo_bits/fn_shift.hpp

@@ -21,6 +21,9 @@
 //! @{
 
 
+// TODO: deprecate shift() in favour of circshift()
+
+
 template<typename T1>
 arma_warn_unused
 arma_inline

include/armadillo_bits/fn_size.hpp

@@ -48,7 +48,7 @@ size(const Base<typename T1::elem_type,T1>& X)
 
 
 
-// explicit overload to workround ADL issues with C++17 std::size()
+// explicit overload to workaround ADL issues with C++17 std::size()
 template<typename eT>
 arma_warn_unused
 inline
@@ -62,7 +62,7 @@ size(const Mat<eT>& X)
 
 
 
-// explicit overload to workround ADL issues with C++17 std::size()
+// explicit overload to workaround ADL issues with C++17 std::size()
 template<typename eT>
 arma_warn_unused
 inline
@@ -76,7 +76,7 @@ size(const Row<eT>& X)
 
 
 
-// explicit overload to workround ADL issues with C++17 std::size()
+// explicit overload to workaround ADL issues with C++17 std::size()
 template<typename eT>
 arma_warn_unused
 inline
@@ -164,7 +164,7 @@ size(const BaseCube<typename T1::elem_type, T1>& X)
 
 
 
-// explicit overload to workround ADL issues with C++17 std::size()
+// explicit overload to workaround ADL issues with C++17 std::size()
 template<typename eT>
 arma_warn_unused
 inline
@@ -242,7 +242,7 @@ size(const SpBase<typename T1::elem_type,T1>& X)
 
 
 
-// explicit overload to workround ADL issues with C++17 std::size()
+// explicit overload to workaround ADL issues with C++17 std::size()
 template<typename eT>
 arma_warn_unused
 inline

include/armadillo_bits/fn_sort_index.hpp

@@ -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));
   }
 
 

include/armadillo_bits/fn_sylvester.hpp

@@ -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);
   }
 
 

include/armadillo_bits/glue_mvnrnd_meat.hpp

@@ -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); } }

include/armadillo_bits/glue_solve_meat.hpp

@@ -180,6 +180,13 @@ glue_solve_gen_full::apply(Mat<eT>& actual_out, const Base<eT,T1>& A_expr, const
     const bool is_sym    = arma_config::optimise_sym && ( (refine || equilibrate || likely_sympd || force_sym || is_band || is_triu || is_tril || auxlib::crippled_lapack(A)) ? false : is_sym_expr<T1>::eval(A_expr.get_ref()) );
     const bool try_sympd = arma_config::optimise_sym && ( (          no_sympd    || is_sym       || force_sym || is_band || is_triu || is_tril || auxlib::crippled_lapack(A)) ? false : (likely_sympd ? true : sym_helper::guess_sympd(A, uword(16))) );
     
+    arma_debug_print("glue_solve_gen_full::apply(): internal flags:");
+    arma_debug_print("is_band:   ", is_band  );
+    arma_debug_print("is_triu:   ", is_triu  );
+    arma_debug_print("is_tril:   ", is_tril  );
+    arma_debug_print("is_sym:    ", is_sym   );
+    arma_debug_print("try_sympd: ", try_sympd);
+    
     if(fast)
       {
       // fast mode: solvers without refinement and without rcond estimate

include/armadillo_bits/glue_times_meat.hpp

@@ -31,20 +31,20 @@ glue_times_redirect2_helper<do_inv_detect>::apply(Mat<typename T1::elem_type>& o
   
   typedef typename T1::elem_type eT;
   
-  const partial_unwrap<T1> tmp1(X.A);
-  const partial_unwrap<T2> tmp2(X.B);
+  const partial_unwrap<T1> U1(X.A);
+  const partial_unwrap<T2> U2(X.B);
   
-  const typename partial_unwrap<T1>::stored_type& A = tmp1.M;
-  const typename partial_unwrap<T2>::stored_type& B = tmp2.M;
+  const typename partial_unwrap<T1>::stored_type& A = U1.M;
+  const typename partial_unwrap<T2>::stored_type& B = U2.M;
   
   constexpr bool use_alpha = partial_unwrap<T1>::do_times || partial_unwrap<T2>::do_times;
-  const     eT       alpha = use_alpha ? (tmp1.get_val() * tmp2.get_val()) : eT(0);
+  const     eT       alpha = use_alpha ? (U1.get_val() * U2.get_val()) : eT(0);
   
   if( (is_cx<eT>::no) && (resolves_to_rowvector<T1>::value && resolves_to_colvector<T2>::value) )
     {
     arma_debug_print("glue_times: dot product optimisation");
     
-    arma_conform_assert_mul_size(A, B, tmp1.do_trans, tmp2.do_trans, "matrix multiplication");
+    arma_conform_assert_mul_size(A, B, U1.do_trans, U2.do_trans, "matrix multiplication");
     
     const eT val = op_dot::direct_dot(A.n_elem, A.memptr(), B.memptr());
     
@@ -55,7 +55,7 @@ glue_times_redirect2_helper<do_inv_detect>::apply(Mat<typename T1::elem_type>& o
     return;
     }
   
-  const bool alias = tmp1.is_alias(out) || tmp2.is_alias(out);
+  const bool alias = U1.is_alias(out) || U2.is_alias(out);
   
   if(alias == false)
     {
@@ -64,7 +64,7 @@ glue_times_redirect2_helper<do_inv_detect>::apply(Mat<typename T1::elem_type>& o
       eT,
       partial_unwrap<T1>::do_trans,
       partial_unwrap<T2>::do_trans,
-      (partial_unwrap<T1>::do_times || partial_unwrap<T2>::do_times)
+      use_alpha
       >
       (out, A, B, alpha);
     }
@@ -77,7 +77,7 @@ glue_times_redirect2_helper<do_inv_detect>::apply(Mat<typename T1::elem_type>& o
       eT,
       partial_unwrap<T1>::do_trans,
       partial_unwrap<T2>::do_trans,
-      (partial_unwrap<T1>::do_times || partial_unwrap<T2>::do_times)
+      use_alpha
       >
       (tmp, A, B, alpha);
     
@@ -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" );
@@ -186,18 +220,18 @@ glue_times_redirect3_helper<do_inv_detect>::apply(Mat<typename T1::elem_type>& o
   // we have exactly 3 objects
   // hence we can safely expand X as X.A.A, X.A.B and X.B
   
-  const partial_unwrap<T1> tmp1(X.A.A);
-  const partial_unwrap<T2> tmp2(X.A.B);
-  const partial_unwrap<T3> tmp3(X.B  );
+  const partial_unwrap<T1> U1(X.A.A);
+  const partial_unwrap<T2> U2(X.A.B);
+  const partial_unwrap<T3> U3(X.B  );
   
-  const typename partial_unwrap<T1>::stored_type& A = tmp1.M;
-  const typename partial_unwrap<T2>::stored_type& B = tmp2.M;
-  const typename partial_unwrap<T3>::stored_type& C = tmp3.M;
+  const typename partial_unwrap<T1>::stored_type& A = U1.M;
+  const typename partial_unwrap<T2>::stored_type& B = U2.M;
+  const typename partial_unwrap<T3>::stored_type& C = U3.M;
   
   constexpr bool use_alpha = partial_unwrap<T1>::do_times || partial_unwrap<T2>::do_times || partial_unwrap<T3>::do_times;
-  const     eT       alpha = use_alpha ? (tmp1.get_val() * tmp2.get_val() * tmp3.get_val()) : eT(0);
+  const     eT       alpha = use_alpha ? (U1.get_val() * U2.get_val() * U3.get_val()) : eT(0);
   
-  const bool alias = tmp1.is_alias(out) || tmp2.is_alias(out) || tmp3.is_alias(out);
+  const bool alias = U1.is_alias(out) || U2.is_alias(out) || U3.is_alias(out);
   
   if(alias == false)
     {
@@ -207,7 +241,7 @@ glue_times_redirect3_helper<do_inv_detect>::apply(Mat<typename T1::elem_type>& o
       partial_unwrap<T1>::do_trans,
       partial_unwrap<T2>::do_trans,
       partial_unwrap<T3>::do_trans,
-      (partial_unwrap<T1>::do_times || partial_unwrap<T2>::do_times || partial_unwrap<T3>::do_times)
+      use_alpha
       >
       (out, A, B, C, alpha);
     }
@@ -221,7 +255,7 @@ glue_times_redirect3_helper<do_inv_detect>::apply(Mat<typename T1::elem_type>& o
       partial_unwrap<T1>::do_trans,
       partial_unwrap<T2>::do_trans,
       partial_unwrap<T3>::do_trans,
-      (partial_unwrap<T1>::do_times || partial_unwrap<T2>::do_times || partial_unwrap<T3>::do_times)
+      use_alpha
       >
       (tmp, A, B, C, alpha);
     
@@ -252,14 +286,14 @@ glue_times_redirect3_helper<true>::apply(Mat<typename T1::elem_type>& out, const
     
     arma_conform_check( (A.is_square() == false), "inv(): given matrix must be square sized" );
     
-    const partial_unwrap<T2> tmp2(X.A.B);
-    const partial_unwrap<T3> tmp3(X.B  );
+    const partial_unwrap<T2> U2(X.A.B);
+    const partial_unwrap<T3> U3(X.B  );
     
-    const typename partial_unwrap<T2>::stored_type& B = tmp2.M;
-    const typename partial_unwrap<T3>::stored_type& C = tmp3.M;
+    const typename partial_unwrap<T2>::stored_type& B = U2.M;
+    const typename partial_unwrap<T3>::stored_type& C = U3.M;
     
     constexpr bool use_alpha = partial_unwrap<T2>::do_times || partial_unwrap<T3>::do_times;
-    const     eT       alpha = use_alpha ? (tmp2.get_val() * tmp3.get_val()) : eT(0);
+    const     eT       alpha = use_alpha ? (U2.get_val() * U3.get_val()) : eT(0);
     
     Mat<eT> BC;
     
@@ -268,7 +302,7 @@ glue_times_redirect3_helper<true>::apply(Mat<typename T1::elem_type>& out, const
       eT,
       partial_unwrap<T2>::do_trans,
       partial_unwrap<T3>::do_trans,
-      (partial_unwrap<T2>::do_times || partial_unwrap<T3>::do_times)
+      use_alpha
       >
       (BC, B, C, alpha);
     
@@ -306,8 +340,8 @@ glue_times_redirect3_helper<true>::apply(Mat<typename T1::elem_type>& out, const
     
     arma_conform_check( (B.is_square() == false), "inv(): given matrix must be square sized" );
     
-    const unwrap<T3> C_tmp(X.B);
-    const Mat<eT>& C = C_tmp.M;
+    const quasi_unwrap<T3> U3(X.B);
+    const Mat<eT>& C =     U3.M;
     
     arma_conform_assert_mul_size(B, C, "matrix multiplication");
     
@@ -330,21 +364,25 @@ glue_times_redirect3_helper<true>::apply(Mat<typename T1::elem_type>& out, const
       return;
       }
     
-    const partial_unwrap_check<T1> tmp1(X.A.A, out);
+    const partial_unwrap<T1> U1(X.A.A);
     
-    const typename partial_unwrap_check<T1>::stored_type& A = tmp1.M;
+    const typename partial_unwrap<T1>::stored_type& A = U1.M;
     
-    constexpr bool use_alpha = partial_unwrap_check<T1>::do_times;
-    const     eT       alpha = use_alpha ? tmp1.get_val() : eT(0);
+    constexpr bool use_alpha = partial_unwrap<T1>::do_times;
+    const     eT       alpha = use_alpha ? U1.get_val() : eT(0);
     
-    glue_times::apply
-      <
-      eT,
-      partial_unwrap_check<T1>::do_trans,
-      false,
-      partial_unwrap_check<T1>::do_times
-      >
-      (out, A, solve_result, alpha);
+    if(U1.is_alias(out))
+      {
+      Mat<eT> tmp;
+      
+      glue_times::apply<eT, partial_unwrap<T1>::do_trans, false, use_alpha>(tmp, A, solve_result, alpha);
+      
+      out.steal_mem(tmp);
+      }
+    else
+      {
+      glue_times::apply<eT, partial_unwrap<T1>::do_trans, false, use_alpha>(out, A, solve_result, alpha);
+      }
     
     return;
     }
@@ -365,16 +403,16 @@ glue_times_redirect<N>::apply(Mat<typename T1::elem_type>& out, const Glue<T1,T2
   
   typedef typename T1::elem_type eT;
   
-  const partial_unwrap<T1> tmp1(X.A);
-  const partial_unwrap<T2> tmp2(X.B);
+  const partial_unwrap<T1> U1(X.A);
+  const partial_unwrap<T2> U2(X.B);
   
-  const typename partial_unwrap<T1>::stored_type& A = tmp1.M;
-  const typename partial_unwrap<T2>::stored_type& B = tmp2.M;
+  const typename partial_unwrap<T1>::stored_type& A = U1.M;
+  const typename partial_unwrap<T2>::stored_type& B = U2.M;
   
   constexpr bool use_alpha = partial_unwrap<T1>::do_times || partial_unwrap<T2>::do_times;
-  const     eT       alpha = use_alpha ? (tmp1.get_val() * tmp2.get_val()) : eT(0);
+  const     eT       alpha = use_alpha ? (U1.get_val() * U2.get_val()) : eT(0);
   
-  const bool alias = tmp1.is_alias(out) || tmp2.is_alias(out);
+  const bool alias = U1.is_alias(out) || U2.is_alias(out);
   
   if(alias == false)
     {
@@ -383,7 +421,7 @@ glue_times_redirect<N>::apply(Mat<typename T1::elem_type>& out, const Glue<T1,T2
       eT,
       partial_unwrap<T1>::do_trans,
       partial_unwrap<T2>::do_trans,
-      (partial_unwrap<T1>::do_times || partial_unwrap<T2>::do_times)
+      use_alpha
       >
       (out, A, B, alpha);
     }
@@ -396,7 +434,7 @@ glue_times_redirect<N>::apply(Mat<typename T1::elem_type>& out, const Glue<T1,T2
       eT,
       partial_unwrap<T1>::do_trans,
       partial_unwrap<T2>::do_trans,
-      (partial_unwrap<T1>::do_times || partial_unwrap<T2>::do_times)
+      use_alpha
       >
       (tmp, A, B, alpha);
     
@@ -446,20 +484,20 @@ glue_times_redirect<4>::apply(Mat<typename T1::elem_type>& out, const Glue< Glue
   // there is exactly 4 objects
   // hence we can safely expand X as X.A.A.A, X.A.A.B, X.A.B and X.B
   
-  const partial_unwrap<T1> tmp1(X.A.A.A);
-  const partial_unwrap<T2> tmp2(X.A.A.B);
-  const partial_unwrap<T3> tmp3(X.A.B  );
-  const partial_unwrap<T4> tmp4(X.B    );
+  const partial_unwrap<T1> U1(X.A.A.A);
+  const partial_unwrap<T2> U2(X.A.A.B);
+  const partial_unwrap<T3> U3(X.A.B  );
+  const partial_unwrap<T4> U4(X.B    );
   
-  const typename partial_unwrap<T1>::stored_type& A = tmp1.M;
-  const typename partial_unwrap<T2>::stored_type& B = tmp2.M;
-  const typename partial_unwrap<T3>::stored_type& C = tmp3.M;
-  const typename partial_unwrap<T4>::stored_type& D = tmp4.M;
+  const typename partial_unwrap<T1>::stored_type& A = U1.M;
+  const typename partial_unwrap<T2>::stored_type& B = U2.M;
+  const typename partial_unwrap<T3>::stored_type& C = U3.M;
+  const typename partial_unwrap<T4>::stored_type& D = U4.M;
   
   constexpr bool use_alpha = partial_unwrap<T1>::do_times || partial_unwrap<T2>::do_times || partial_unwrap<T3>::do_times || partial_unwrap<T4>::do_times;
-  const     eT       alpha = use_alpha ? (tmp1.get_val() * tmp2.get_val() * tmp3.get_val() * tmp4.get_val()) : eT(0);
+  const     eT       alpha = use_alpha ? (U1.get_val() * U2.get_val() * U3.get_val() * U4.get_val()) : eT(0);
   
-  const bool alias = tmp1.is_alias(out) || tmp2.is_alias(out) || tmp3.is_alias(out) || tmp4.is_alias(out);
+  const bool alias = U1.is_alias(out) || U2.is_alias(out) || U3.is_alias(out) || U4.is_alias(out);
   
   if(alias == false)
     {
@@ -470,7 +508,7 @@ glue_times_redirect<4>::apply(Mat<typename T1::elem_type>& out, const Glue< Glue
       partial_unwrap<T2>::do_trans,
       partial_unwrap<T3>::do_trans,
       partial_unwrap<T4>::do_trans,
-      (partial_unwrap<T1>::do_times || partial_unwrap<T2>::do_times || partial_unwrap<T3>::do_times || partial_unwrap<T4>::do_times)
+      use_alpha
       >
       (out, A, B, C, D, alpha);
     }
@@ -485,7 +523,7 @@ glue_times_redirect<4>::apply(Mat<typename T1::elem_type>& out, const Glue< Glue
       partial_unwrap<T2>::do_trans,
       partial_unwrap<T3>::do_trans,
       partial_unwrap<T4>::do_trans,
-      (partial_unwrap<T1>::do_times || partial_unwrap<T2>::do_times || partial_unwrap<T3>::do_times || partial_unwrap<T4>::do_times)
+      use_alpha
       >
       (tmp, A, B, C, D, alpha);
     
@@ -504,7 +542,7 @@ glue_times::apply(Mat<typename T1::elem_type>& out, const Glue<T1,T2,glue_times>
   
   constexpr uword N_mat = 1 + depth_lhs< glue_times, Glue<T1,T2,glue_times> >::num;
   
-  arma_debug_print(arma_str::format("N_mat: %u") % N_mat);
+  arma_debug_print(arma_str::format("glue_times::apply(): N_mat: %u") % N_mat);
   
   glue_times_redirect<N_mat>::apply(out, X);
   }
@@ -518,7 +556,11 @@ glue_times::apply_inplace(Mat<typename T1::elem_type>& out, const T1& X)
   {
   arma_debug_sigprint();
   
-  out = out * X;
+  typedef typename T1::elem_type eT;
+  
+  Mat<eT> tmp = out * X;
+  
+  out.steal_mem(tmp);
   }
 
 
@@ -533,9 +575,9 @@ glue_times::apply_inplace_plus(Mat<typename T1::elem_type>& out, const Glue<T1,
   typedef typename T1::elem_type            eT;
   typedef typename get_pod_type<eT>::result  T;
   
-  if( (is_outer_product<T1>::value) || (has_op_inv_any<T1>::value) || (has_op_inv_any<T2>::value) )
+  if( X.is_alias(out) || (is_outer_product<T1>::value) || (has_op_inv_any<T1>::value) || (has_op_inv_any<T2>::value) )
     {
-    // partial workaround for corner cases
+    // handle aliasing and partial workaround for corner cases
     
     const Mat<eT> tmp(X);
     
@@ -544,21 +586,21 @@ glue_times::apply_inplace_plus(Mat<typename T1::elem_type>& out, const Glue<T1,
     return;
     }
   
-  const partial_unwrap_check<T1> tmp1(X.A, out);
-  const partial_unwrap_check<T2> tmp2(X.B, out);
+  const partial_unwrap<T1> U1(X.A);
+  const partial_unwrap<T2> U2(X.B);
   
-  typedef typename partial_unwrap_check<T1>::stored_type TA;
-  typedef typename partial_unwrap_check<T2>::stored_type TB;
+  typedef typename partial_unwrap<T1>::stored_type TA;
+  typedef typename partial_unwrap<T2>::stored_type TB;
   
-  const TA& A = tmp1.M;
-  const TB& B = tmp2.M;
+  const TA& A = U1.M;
+  const TB& B = U2.M;
   
-  constexpr bool do_trans_A = partial_unwrap_check<T1>::do_trans;
-  constexpr bool do_trans_B = partial_unwrap_check<T2>::do_trans;
+  constexpr bool do_trans_A = partial_unwrap<T1>::do_trans;
+  constexpr bool do_trans_B = partial_unwrap<T2>::do_trans;
   
-  const bool use_alpha = partial_unwrap_check<T1>::do_times || partial_unwrap_check<T2>::do_times || (sign < sword(0));
+  const bool use_alpha = partial_unwrap<T1>::do_times || partial_unwrap<T2>::do_times || (sign < sword(0));
   
-  const eT       alpha = use_alpha ? ( tmp1.get_val() * tmp2.get_val() * ( (sign > sword(0)) ? eT(1) : eT(-1) ) ) : eT(0);
+  const eT       alpha = use_alpha ? ( U1.get_val() * U2.get_val() * ( (sign > sword(0)) ? eT(1) : eT(-1) ) ) : eT(0);
   
   arma_conform_assert_mul_size(A, B, do_trans_A, do_trans_B, "matrix multiplication");
   

include/armadillo_bits/gmm_diag_meat.hpp

@@ -2239,7 +2239,7 @@ gmm_diag<eT>::km_iterate(const Mat<eT>& X, const uword max_iter, const bool verb
           }
         else
           {
-          // recover by using a randomly seleced sample (last resort)
+          // recover by using a randomly selected sample (last resort)
           proposed_i = as_scalar(randi<uvec>(1, distr_param(0,X_n_cols-1)));
           }
         
@@ -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; }
       }
     
     

include/armadillo_bits/gmm_full_meat.hpp

@@ -937,7 +937,7 @@ gmm_full<eT>::init_constants(const bool calc_chol)
           Mat<eT>& inv_fcov =  inv_fcovs.slice(g);
     
   //const bool inv_ok = auxlib::inv(tmp_inv, fcov);
-    const bool inv_ok = auxlib::inv_sympd(tmp_inv, fcov);
+    const bool inv_ok = auxlib::inv_sympd(tmp_inv, fcov);  // using inv_sympd() instead of inv() to ensure we can do cholesky decomp via op_chol
     
     eT log_det_val  = eT(0);
     eT log_det_sign = eT(0);
@@ -2270,7 +2270,7 @@ gmm_full<eT>::km_iterate(const Mat<eT>& X, const uword max_iter, const bool verb
           }
         else
           {
-          // recover by using a randomly seleced sample (last resort)
+          // recover by using a randomly selected sample (last resort)
           proposed_i = as_scalar(randi<uvec>(1, distr_param(0,X_n_cols-1)));
           }
         
@@ -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);
     

include/armadillo_bits/include_superlu.hpp

@@ -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
@@ -56,7 +56,7 @@
 
 #undef ARMA_SLU_HEADERS_FOUND
 
-// Since we need to suport float, double, cx_float and cx_double,
+// Since we need to support float, double, cx_float and cx_double,
 // as well as preserve the sanity of the user,
 // we cannot simply include all the SuperLU headers due to their messy state
 // (duplicate definitions, pollution of global namespace, bizarro defines).
@@ -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;

include/armadillo_bits/memory.hpp

@@ -90,7 +90,7 @@ memory::acquire(const uword n_elem)
   #else
     {
     //return ( new(std::nothrow) eT[n_elem] );
-    out_memptr = (eT *) malloc(sizeof(eT)*n_elem);
+    out_memptr = (eT *) std::malloc(sizeof(eT)*n_elem);
     }
   #endif
   
@@ -124,7 +124,7 @@ memory::release(eT* mem)
     }
   #elif defined(ARMA_HAVE_POSIX_MEMALIGN)
     {
-    free( (void *)(mem) );
+    std::free( (void *)(mem) );
     }
   #elif defined(_MSC_VER)
     {
@@ -134,7 +134,7 @@ memory::release(eT* mem)
   #else
     {
     //delete [] mem;
-    free( (void *)(mem) );
+    std::free( (void *)(mem) );
     }
   #endif
   

include/armadillo_bits/mtGlueCube_bones.hpp

@@ -33,9 +33,9 @@ class mtGlueCube : public BaseCube< out_eT, mtGlueCube<out_eT, T1, T2, glue_type
   arma_inline  mtGlueCube(const T1& in_A, const T2& in_B, const uword in_aux_uword);
   arma_inline ~mtGlueCube();
   
-  arma_aligned const T1&   A;         //!< first operand;  must be derived from BaseCube
-  arma_aligned const T2&   B;         //!< second operand; must be derived from BaseCube
-  arma_aligned       uword aux_uword; //!< storage of auxiliary data, uword format
+  const T1&   A;         //!< first operand;  must be derived from BaseCube
+  const T2&   B;         //!< second operand; must be derived from BaseCube
+        uword aux_uword; //!< storage of auxiliary data, uword format
   };
 
 

include/armadillo_bits/mtGlue_bones.hpp

@@ -40,9 +40,9 @@ class mtGlue : public Base< out_eT, mtGlue<out_eT, T1, T2, glue_type> >
   template<typename eT2>
   inline bool is_alias(const Mat<eT2>& X) const;
   
-  arma_aligned const T1&   A;         //!< first operand;  must be derived from Base
-  arma_aligned const T2&   B;         //!< second operand; must be derived from Base
-  arma_aligned       uword aux_uword; //!< storage of auxiliary data, uword format
+  const T1&   A;         //!< first operand;  must be derived from Base
+  const T2&   B;         //!< second operand; must be derived from Base
+        uword aux_uword; //!< storage of auxiliary data, uword format
   };