patch bump

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


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


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

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


 # extract version from sources
@@ -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/)
@@ -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)
@@ -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)

 ---

@@ -303,7 +303,7 @@ the active solution platform must be set to x64, instead of win32.

 The MSVC project was tested on Windows 10 (64 bit) with Visual Studio C++ 2019.
 Adaptations may be required for 32 bit systems, later versions of Windows and/or the compiler.
-For example, options such as `ARMA_BLAS_LONG` and `ARMA_BLAS_UNDERSCORE`,
+For example, options such as `ARMA_BLAS_LONG_LONG` and `ARMA_BLAS_UNDERSCORE`,
 defined in `include/armadillo_bits/config.hpp`, may need to be either enabled or disabled.

 The folder `examples/lib_win64` contains a copy of lib and dll files
@@ -333,29 +333,28 @@ as used by the replacement libraries. Specifically, the following defines
 may need to be enabled or disabled:

    ARMA_USE_WRAPPER  
-    ARMA_BLAS_CAPITALS  
-    ARMA_BLAS_UNDERSCORE  
-    ARMA_BLAS_LONG  
    ARMA_BLAS_LONG_LONG  
+    ARMA_DONT_USE_FORTRAN_HIDDEN_ARGS  
+    ARMA_BLAS_UNDERSCORE  
+    ARMA_BLAS_CAPITALS  

 See the documentation for more information on the above defines.

 On Linux-based systems, MKL might be installed in a non-standard location such as `/opt`
 which can cause problems during linking.
-Before installing Armadillo, the system should know where the MKL libraries are located.
-For example, `/opt/intel/mkl/lib/intel64/`.
-This can be achieved by setting the `LD_LIBRARY_PATH` environment variable,
-or for a more permanent solution, adding the directory locations to `/etc/ld.so.conf`.
-It may also be possible to store a text file with the locations
-in the `/etc/ld.so.conf.d` directory. For example, `/etc/ld.so.conf.d/mkl.conf`.
-If `/etc/ld.so.conf` is modified or `/etc/ld.so.conf.d/mkl.conf` is created,
-`/sbin/ldconfig` must be run afterwards.
+Examples: `/opt/intel/oneapi/mkl/latest/lib`, `/opt/intel/mkl/lib/intel64/`.

-Below is an example of `/etc/ld.so.conf.d/mkl.conf`
-where Intel MKL is installed in `/opt/intel`
+Before installing Armadillo, the system must know where the MKL libraries are located.
+This can be achieved via several ways:

-    /opt/intel/lib/intel64  
-    /opt/intel/mkl/lib/intel64  
+1. By setting the `LD_LIBRARY_PATH` environment variable.
+
+2. By adding the MKL library directory locations to the `/etc/ld.so.conf` text file,
+   followed by running `/sbin/ldconfig`.
+
+3. By creating a text file named `/etc/ld.so.conf.d/mkl.conf`
+   which contains the MKL library directory locations,
+   followed by running `/sbin/ldconfig`.

 If MKL is installed and it is persistently giving problems during linking,
 Support for MKL can be disabled by editing the CMakeLists.txt file,
@@ -370,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`.

 ---
@@ -461,12 +460,12 @@ Octave/Matlab with C++ code that uses Armadillo matrices.

 ### 16: Related Software Using Armadillo

-* ensmallen: C++ library for non-linear numerical optimisation (L-BFGS, SGD, CMA-ES, etc)  
-  https://ensmallen.org/
-
 * MLPACK: extensive library of machine learning algorithms  
  https://mlpack.org

+* ensmallen: C++ library for numerical optimisation (L-BFGS, SGD, CMA-ES, etc)  
+  https://ensmallen.org/
+
 * RcppArmadillo: integration of Armadillo with R  
  https://dirk.eddelbuettel.com/code/rcpp.armadillo.html

@@ -115,6 +115,9 @@ foreach (MKL_NAME ${MKL_NAMES})
    /usr/local/mkl/*/lib/em64t
    /usr/local/mkl/*/lib/lib64
    /usr/local/mkl/*/lib/32
+    /usr/local/intel/oneapi/mkl/latest/lib
+    /opt/intel/oneapi/mkl/latest/lib
+    /opt/intel/oneapi/mkl/*/lib
    )

  set(TMP_LIBRARY ${${MKL_NAME}_LIBRARY})
@@ -0,0 +1,97 @@
+# Try to find suitable version of SuperLU
+# Once done, this will define
+#
+#  SuperLU_FOUND        - system has SuperLU
+#  SuperLU_LIBRARY      - Link this to use SuperLU
+#  SuperLU_INCLUDE_DIR  - directory of SuperLU headers
+
+find_path(SuperLU_INCLUDE_DIR slu_ddefs.h
+  PATHS ${CMAKE_SYSTEM_INCLUDE_PATH} /usr/include /usr/local/include /opt/local/include
+  PATH_SUFFIXES superlu SuperLU ""
+)
+
+find_library(SuperLU_LIBRARY
+  NAMES superlu
+  PATHS ${CMAKE_SYSTEM_LIBRARY_PATH} /usr/lib64 /usr/lib /usr/local/lib64 /usr/local/lib /opt/local/lib64 /opt/local/lib
+)
+
+set(SuperLU_FOUND NO)
+
+if(NOT SuperLU_FIND_QUIETLY)
+  message(STATUS "Looking for suitable version of SuperLU")
+endif()
+
+if(SuperLU_LIBRARY AND SuperLU_INCLUDE_DIR)
+  
+  # determine if this is a suitable version of SuperLU
+  
+  if(EXISTS "${SuperLU_INCLUDE_DIR}/slu_ddefs.h")
+    
+    file(READ "${SuperLU_INCLUDE_DIR}/slu_ddefs.h" _slu_HEADER_CONTENTS)
+    
+    string(REGEX REPLACE ".*version ([0-9]).*" "\\1" SLU_VERSION_MAJOR "${_slu_HEADER_CONTENTS}")
+    
+    if("${SLU_VERSION_MAJOR}" EQUAL "7")
+      set(SuperLU_FOUND YES)
+    endif()
+    
+    if("${SLU_VERSION_MAJOR}" EQUAL "6")
+      set(SuperLU_FOUND YES)
+    endif()
+    
+    if("${SLU_VERSION_MAJOR}" EQUAL "5")
+      set(SuperLU_FOUND YES)
+    endif()
+    
+    if("${SLU_VERSION_MAJOR}" EQUAL "4")
+      
+      # NOTE: SLU_VERSION_MAJOR = 4 is not 100% reliable, as the version string in slu_ddefs.h wasn't updated in SuperLU 5.0, 5.1, 5.2
+      # NOTE: We need to explicitly check the prototype of dgssvx() to determine whether this is version 4 or 5.
+      
+      string(FIND "${_slu_HEADER_CONTENTS}" "dgssvx" _slu_dgssvx_start)
+      
+      if("${_slu_dgssvx_start}" GREATER 0)
+        
+        string(SUBSTRING "${_slu_HEADER_CONTENTS}" ${_slu_dgssvx_start} -1 _slu_dgssvx_tmp)
+        
+        string(FIND "${_slu_dgssvx_tmp}" ";" _slu_dgssvx_len)
+        
+        string(SUBSTRING "${_slu_HEADER_CONTENTS}" ${_slu_dgssvx_start} ${_slu_dgssvx_len} _slu_dgssvx_proto)
+        
+        string(FIND "${_slu_dgssvx_proto}" "GlobalLU_t" _slu_dgssvx_check)
+        
+        if("${_slu_dgssvx_check}" EQUAL -1)
+          # in version 4, dgssvx() doesn't have a GlobalLU_t argument
+          if(NOT SuperLU_FIND_QUIETLY)
+            message(STATUS "Found SuperLU 4, which is insufficient")
+          endif()
+        else()
+          set(SuperLU_FOUND YES)
+        endif()
+      
+      endif()
+    
+    endif()
+  
+  else()
+  
+    message(STATUS "Could not read ${SuperLU_INCLUDE_DIR}/slu_ddefs.h to verify suitable version")
+  
+  endif()
+  
+endif()
+
+if(SuperLU_LIBRARY AND NOT SuperLU_INCLUDE_DIR)
+  if(NOT SuperLU_FIND_QUIETLY)
+    message(STATUS "Incomplete installation of SuperLU: no SuperLU headers found")
+  endif()
+endif()
+  
+if(NOT SuperLU_FOUND) 
+  if(NOT SuperLU_FIND_QUIETLY)
+    message(STATUS "Could not find suitable version of SuperLU")
+  endif()
+else()
+  message(STATUS "Found SuperLU library: ${SuperLU_LIBRARY}")
+endif()
+
@@ -1,90 +0,0 @@
-# Try to find SuperLU 5.x
-# Once done, this will define
-#
-#  SuperLU_FOUND        - system has SuperLU
-#  SuperLU_LIBRARY      - Link this to use SuperLU
-#  SuperLU_INCLUDE_DIR  - directory of SuperLU headers
-
-find_path(SuperLU_INCLUDE_DIR slu_ddefs.h
-  PATHS ${CMAKE_SYSTEM_INCLUDE_PATH} /usr/include /usr/local/include /opt/local/include
-  PATH_SUFFIXES superlu SuperLU ""
-)
-
-find_library(SuperLU_LIBRARY
-  NAMES superlu
-  PATHS ${CMAKE_SYSTEM_LIBRARY_PATH} /usr/lib64 /usr/lib /usr/local/lib64 /usr/local/lib /opt/local/lib64 /opt/local/lib
-)
-
-set(SuperLU_FOUND NO)
-
-if(NOT SuperLU_FIND_QUIETLY)
-  message(STATUS "Looking for SuperLU version 5")
-endif()
-
-if(SuperLU_LIBRARY AND SuperLU_INCLUDE_DIR)
-  
-  # Now determine if the version is version 5.
-  # This voodoo is required as SuperLU developers like to keep things messy.
-  
-  if(EXISTS "${SuperLU_INCLUDE_DIR}/slu_ddefs.h")
-    
-    file(READ "${SuperLU_INCLUDE_DIR}/slu_ddefs.h" _slu_HEADER_CONTENTS)
-    
-    string(REGEX REPLACE ".*version ([0-9]).*" "\\1" SLU_VERSION_MAJOR "${_slu_HEADER_CONTENTS}")
-    
-    if("${SLU_VERSION_MAJOR}" EQUAL "5")
-      set(SuperLU_FOUND YES)
-    else()
-      if("${SLU_VERSION_MAJOR}" EQUAL "4")
-        
-        # NOTE: SLU_VERSION_MAJOR = 4 is not 100% reliable, as the version string in slu_ddefs.h wasn't updated in SuperLU 5.0, 5.1, 5.2
-        # NOTE: We need to explicitly check the prototype of dgssvx() to determine whether this is version 4 or 5.
-        
-        string(FIND "${_slu_HEADER_CONTENTS}" "dgssvx" _slu_dgssvx_start)
-        
-        if("${_slu_dgssvx_start}" GREATER 0)
-          
-          string(SUBSTRING "${_slu_HEADER_CONTENTS}" ${_slu_dgssvx_start} -1 _slu_dgssvx_tmp)
-          
-          string(FIND "${_slu_dgssvx_tmp}" ";" _slu_dgssvx_len)
-          
-          string(SUBSTRING "${_slu_HEADER_CONTENTS}" ${_slu_dgssvx_start} ${_slu_dgssvx_len} _slu_dgssvx_proto)
-          
-          string(FIND "${_slu_dgssvx_proto}" "GlobalLU_t" _slu_dgssvx_check)
-          
-          if("${_slu_dgssvx_check}" EQUAL -1)
-            # in version 4, dgssvx() doesn't have a GlobalLU_t argument
-            if(NOT SuperLU_FIND_QUIETLY)
-              message(STATUS "Found SuperLU, but it doesn't appear to be version 5")
-            endif()
-          else()
-            set(SuperLU_FOUND YES)
-          endif()
-        endif()
-      else()
-        if(NOT SuperLU_FIND_QUIETLY)
-          message(STATUS "Found SuperLU version ${SLU_VERSION_MAJOR}, but version 5 is required")
-        endif()
-      endif()
-    endif()
-  else()
-    if(NOT SuperLU_FIND_QUIETLY)
-      message(STATUS "Found SuperLU and headers, but could not verify version 5")
-    endif()
-  endif()
-endif()
-
-if(SuperLU_LIBRARY AND NOT SuperLU_INCLUDE_DIR)
-  if(NOT SuperLU_FIND_QUIETLY)
-    message(STATUS "Incomplete installation of SuperLU: no SuperLU headers found")
-  endif()
-endif()
-  
-if(NOT SuperLU_FOUND) 
-  if(NOT SuperLU_FIND_QUIETLY)
-    message(STATUS "Could not find SuperLU")
-  endif()
-else()
-  message(STATUS "Found SuperLU: ${SuperLU_LIBRARY}")
-endif()
-
@@ -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:
@@ -304,6 +304,7 @@ namespace arma
  #include "armadillo_bits/op_clamp_bones.hpp"
  #include "armadillo_bits/op_expmat_bones.hpp"
  #include "armadillo_bits/op_nonzeros_bones.hpp"
+  #include "armadillo_bits/op_omit_bones.hpp"
  #include "armadillo_bits/op_diff_bones.hpp"
  #include "armadillo_bits/op_norm_bones.hpp"
  #include "armadillo_bits/op_vecnorm_bones.hpp"
@@ -377,6 +378,7 @@ namespace arma
  #include "armadillo_bits/spop_norm_bones.hpp"
  #include "armadillo_bits/spop_shift_bones.hpp"
  #include "armadillo_bits/spop_relational_bones.hpp"
+  #include "armadillo_bits/spop_omit_bones.hpp"
  
  #include "armadillo_bits/spglue_plus_bones.hpp"
  #include "armadillo_bits/spglue_minus_bones.hpp"
@@ -564,6 +566,7 @@ namespace arma
  #include "armadillo_bits/fn_clamp.hpp"
  #include "armadillo_bits/fn_expmat.hpp"
  #include "armadillo_bits/fn_nonzeros.hpp"
+  #include "armadillo_bits/fn_omit.hpp"
  #include "armadillo_bits/fn_interp1.hpp"
  #include "armadillo_bits/fn_interp2.hpp"
  #include "armadillo_bits/fn_qz.hpp"
@@ -590,6 +593,7 @@ namespace arma
  #include "armadillo_bits/fn_powmat.hpp"
  #include "armadillo_bits/fn_powext.hpp"
  #include "armadillo_bits/fn_diags_spdiags.hpp"
+  #include "armadillo_bits/fn_balance.hpp"
  
  #include "armadillo_bits/fn_speye.hpp"
  #include "armadillo_bits/fn_spones.hpp"
@@ -758,6 +762,7 @@ namespace arma
  #include "armadillo_bits/op_clamp_meat.hpp"
  #include "armadillo_bits/op_expmat_meat.hpp"
  #include "armadillo_bits/op_nonzeros_meat.hpp"
+  #include "armadillo_bits/op_omit_meat.hpp"
  #include "armadillo_bits/op_diff_meat.hpp"
  #include "armadillo_bits/op_norm_meat.hpp"
  #include "armadillo_bits/op_vecnorm_meat.hpp"
@@ -831,6 +836,7 @@ namespace arma
  #include "armadillo_bits/spop_norm_meat.hpp"
  #include "armadillo_bits/spop_shift_meat.hpp"
  #include "armadillo_bits/spop_relational_meat.hpp"
+  #include "armadillo_bits/spop_omit_meat.hpp"
  
  #include "armadillo_bits/spglue_plus_meat.hpp"
  #include "armadillo_bits/spglue_minus_meat.hpp"
@@ -328,7 +328,7 @@ BaseCube<elem_type,derived>::is_finite() const
    for(uword c=0; c<n_c; ++c)
    for(uword r=0; r<n_r; ++r)
      {
-      if( arma_isfinite(P.at(r,c,s)) == false )  { return false; }
+      if( arma_isnonfinite(P.at(r,c,s)) )  { return false; }
      }
    }
  
@@ -436,7 +436,7 @@ BaseCube<elem_type,derived>::has_nonfinite() const
    for(uword c=0; c<n_c; ++c)
    for(uword r=0; r<n_r; ++r)
      {
-      if(arma_isfinite(P.at(r,c,s)) == false)  { return true; }
+      if(arma_isnonfinite(P.at(r,c,s)))  { return true; }
      }
    }
  
@@ -127,11 +127,11 @@ struct Base
  arma_warn_unused inline elem_type min() const;
  arma_warn_unused inline elem_type max() const;
  
-  inline elem_type min(uword& index_of_min_val) const;
-  inline elem_type max(uword& index_of_max_val) const;
+  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;
  
-  inline elem_type min(uword& row_of_min_val, uword& col_of_min_val) const;
-  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;
@@ -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; }
        }
      }
    }
@@ -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 ];
@@ -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(); }
  }


@@ -33,8 +33,11 @@ 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; }
  
  static constexpr bool is_row  = op_type::template traits<T1>::is_row;
  static constexpr bool is_col  = op_type::template traits<T1>::is_col;
@@ -387,12 +387,11 @@ class Cube : public BaseCube< eT, Cube<eT> >
  arma_warn_unused inline eT min() const;
  arma_warn_unused inline eT max() const;
  
-  inline eT min(uword& index_of_min_val) const;
-  inline eT max(uword& index_of_max_val) const;
-  
-  inline eT min(uword& row_of_min_val, uword& col_of_min_val, uword& slice_of_min_val) const;
-  inline eT max(uword& row_of_max_val, uword& col_of_max_val, uword& slice_of_max_val) const;
+  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_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;
@@ -486,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;
  
@@ -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" );
  }
@@ -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();
@@ -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();  }
@@ -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();
@@ -54,6 +54,9 @@ class Gen
  inline void apply_inplace_div  (Mat<elem_type>& out) const;
  
  inline void apply(subview<elem_type>& out) const;
+  
+  template<typename eT2>
+  constexpr bool is_alias(const Mat<eT2>&) const { return false; }
  };


@@ -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();  }
@@ -56,6 +56,9 @@ class Glue
  inline  Glue(const T1& in_A, const T2& in_B, const uword in_aux_uword);
  inline ~Glue();
  
+  template<typename eT2>
+  inline bool is_alias(const Mat<eT2>& X) const;
+  
  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
@@ -53,4 +53,17 @@ Glue<T1,T2,glue_type>::~Glue()



+template<typename T1, typename T2, typename glue_type>
+template<typename eT2>
+inline
+bool
+Glue<T1,T2,glue_type>::is_alias(const Mat<eT2>& X) const
+  {
+  arma_debug_sigprint();
+  
+  return (A.is_alias(X) || B.is_alias(X));
+  }
+
+
+
 //! @}
@@ -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();
        }
      }
    }
@@ -541,11 +541,11 @@ class Mat : public Base< eT, Mat<eT> >
  arma_warn_unused inline eT min() const;
  arma_warn_unused inline eT max() const;
  
-  inline eT min(uword& index_of_min_val) const;
-  inline eT max(uword& index_of_max_val) const;
+  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;
  
-  inline eT min(uword& row_of_min_val, uword& col_of_min_val) const;
-  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;
@@ -771,6 +771,9 @@ class Mat : public Base< eT, Mat<eT> >
  
  inline void steal_mem_col(Mat& X, const uword max_n_rows);
  
+  template<typename eT2>
+  arma_inline bool is_alias(const Mat<eT2>& X) const; //!< don't use this unless you're writing code internal to Armadillo
+  
  
  template<uword fixed_n_rows, uword fixed_n_cols> class fixed;
  
@@ -822,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);
  
@@ -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;
  }
@@ -1324,6 +1328,19 @@ Mat<eT>::steal_mem_col(Mat<eT>& x, const uword max_n_rows)



+template<typename eT>
+template<typename eT2>
+arma_inline
+bool
+Mat<eT>::is_alias(const Mat<eT2>& X) const
+  {
+  arma_debug_sigprint();
+  
+  return (is_same_type<eT,eT2>::yes) && (void_ptr(this) == void_ptr(&X));
+  }
+
+
+
 //! construct a matrix from a given auxiliary array of eTs.
 //! if copy_aux_mem is true, new memory is allocated and the array is copied.
 //! if copy_aux_mem is false, the auxiliary array is used directly (without allocating memory and copying).
@@ -2054,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>&
@@ -2069,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>&
@@ -2244,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>&
@@ -2261,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>&
@@ -2370,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>&
@@ -2385,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>&
@@ -2849,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;
@@ -2876,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();
    
@@ -2900,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;
@@ -2927,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();
    
@@ -2951,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;
@@ -2978,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();
    
@@ -5172,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);
  }

@@ -5194,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;
@@ -5215,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;
@@ -5236,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;
@@ -5274,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;
@@ -5295,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;
@@ -6378,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&
@@ -6391,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&
@@ -6447,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&
@@ -6460,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&
@@ -6613,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
@@ -9552,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(); }
  }


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


@@ -58,10 +58,13 @@ class Op
  inline          Op(const T1& in_m, const uword     in_aux_uword_a, const uword in_aux_uword_b);
  inline         ~Op();
  
-  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
+  template<typename eT2>
+  inline bool is_alias(const Mat<eT2>& X) const;
+  
+  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
  };


@@ -76,4 +76,17 @@ Op<T1, op_type>::~Op()



+template<typename T1, typename op_type>
+template<typename eT2>
+inline
+bool
+Op<T1, op_type>::is_alias(const Mat<eT2>& X) const
+  {
+  arma_debug_sigprint();
+  
+  return m.is_alias(X);
+  }
+
+
+
 //! @}
@@ -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)
@@ -188,7 +188,7 @@ struct Proxy< Mat<eT> >
  arma_inline aligned_ea_type get_aligned_ea() const { return Q;          }
  
  template<typename eT2>
-  arma_inline bool is_alias(const Mat<eT2>& X) const { return (is_same_type<eT,eT2>::value) ? (void_ptr(&Q) == void_ptr(&X)) : false; }
+  arma_inline bool is_alias(const Mat<eT2>& X) const { return (is_same_type<eT,eT2>::yes) && (void_ptr(&Q) == void_ptr(&X)); }
  
  template<typename eT2>
  arma_inline bool has_overlap(const subview<eT2>& X) const { return is_alias(X.m); }
@@ -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)
@@ -235,7 +235,7 @@ struct Proxy< Col<eT> >
  arma_inline aligned_ea_type get_aligned_ea() const { return Q;          }
  
  template<typename eT2>
-  arma_inline bool is_alias(const Mat<eT2>& X) const { return (is_same_type<eT,eT2>::value) ? (void_ptr(&Q) == void_ptr(&X)) : false; }
+  arma_inline bool is_alias(const Mat<eT2>& X) const { return (is_same_type<eT,eT2>::yes) && (void_ptr(&Q) == void_ptr(&X)); }
  
  template<typename eT2>
  arma_inline bool has_overlap(const subview<eT2>& X) const { return is_alias(X.m); }
@@ -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)
@@ -282,7 +282,7 @@ struct Proxy< Row<eT> >
  arma_inline aligned_ea_type get_aligned_ea() const { return Q;          }
  
  template<typename eT2>
-  arma_inline bool is_alias(const Mat<eT2>& X) const { return (is_same_type<eT,eT2>::value) ? (void_ptr(&Q) == void_ptr(&X)) : false; }
+  arma_inline bool is_alias(const Mat<eT2>& X) const { return (is_same_type<eT,eT2>::yes) && (void_ptr(&Q) == void_ptr(&X)); }
  
  template<typename eT2>
  arma_inline bool has_overlap(const subview<eT2>& X) const { return is_alias(X.m); }
@@ -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)
@@ -1013,7 +1013,7 @@ struct Proxy< subview<eT> >
  arma_inline aligned_ea_type get_aligned_ea() const { return Q; }
  
  template<typename eT2>
-  arma_inline bool is_alias(const Mat<eT2>& X) const { return (is_same_type<eT,eT2>::value) ? (void_ptr(&(Q.m)) == void_ptr(&X)) : false; }
+  arma_inline bool is_alias(const Mat<eT2>& X) const { return (is_same_type<eT,eT2>::yes) && (void_ptr(&(Q.m)) == void_ptr(&X)); }
  
  template<typename eT2>
  arma_inline bool has_overlap(const subview<eT2>& X) const { return Q.check_overlap(X); }
@@ -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)
@@ -1060,7 +1060,7 @@ struct Proxy< subview_col<eT> >
  arma_inline aligned_ea_type get_aligned_ea() const { return Q;        }
  
  template<typename eT2>
-  arma_inline bool is_alias(const Mat<eT2>& X) const { return (is_same_type<eT,eT2>::value) ? (void_ptr(&(Q.m)) == void_ptr(&X)) : false; }
+  arma_inline bool is_alias(const Mat<eT2>& X) const { return (is_same_type<eT,eT2>::yes) && (void_ptr(&(Q.m)) == void_ptr(&X)); }
  
  template<typename eT2>
  arma_inline bool has_overlap(const subview<eT2>& X) const { return Q.check_overlap(X); }
@@ -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)
@@ -1109,7 +1109,7 @@ struct Proxy< subview_cols<eT> >
  arma_inline aligned_ea_type get_aligned_ea() const { return Q;          }
  
  template<typename eT2>
-  arma_inline bool is_alias(const Mat<eT2>& X) const { return (is_same_type<eT,eT2>::value) ? (void_ptr(&(sv.m)) == void_ptr(&X)) : false; }
+  arma_inline bool is_alias(const Mat<eT2>& X) const { return (is_same_type<eT,eT2>::yes) && (void_ptr(&(sv.m)) == void_ptr(&X)); }
  
  template<typename eT2>
  arma_inline bool has_overlap(const subview<eT2>& X) const { return sv.check_overlap(X); }
@@ -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)
@@ -1156,7 +1156,7 @@ struct Proxy< subview_row<eT> >
  arma_inline aligned_ea_type get_aligned_ea() const { return Q; }
  
  template<typename eT2>
-  arma_inline bool is_alias(const Mat<eT2>& X) const { return (is_same_type<eT,eT2>::value) ? (void_ptr(&(Q.m)) == void_ptr(&X)) : false; }
+  arma_inline bool is_alias(const Mat<eT2>& X) const { return (is_same_type<eT,eT2>::yes) && (void_ptr(&(Q.m)) == void_ptr(&X)); }
  
  template<typename eT2>
  arma_inline bool has_overlap(const subview<eT2>& X) const { return Q.check_overlap(X); }
@@ -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)
@@ -1304,7 +1304,7 @@ struct Proxy< diagview<eT> >
  arma_inline aligned_ea_type get_aligned_ea() const { return Q; }
  
  template<typename eT2>
-  arma_inline bool is_alias(const Mat<eT2>& X) const { return (is_same_type<eT,eT2>::value) ? (void_ptr(&(Q.m)) == void_ptr(&X)) : false; }
+  arma_inline bool is_alias(const Mat<eT2>& X) const { return (is_same_type<eT,eT2>::yes) && (void_ptr(&(Q.m)) == void_ptr(&X)); }
  
  template<typename eT2>
  arma_inline bool has_overlap(const subview<eT2>& X) const { return is_alias(X.m); }
@@ -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)
@@ -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)
@@ -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 ];
@@ -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(); }
  }


@@ -76,11 +76,11 @@ struct SpBase
  arma_warn_unused inline elem_type min() const;
  arma_warn_unused inline elem_type max() const;
  
-  inline elem_type min(uword& index_of_min_val) const;
-  inline elem_type max(uword& index_of_max_val) const;
+  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;
  
-  inline elem_type min(uword& row_of_min_val, uword& col_of_min_val) const;
-  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;
@@ -748,7 +748,7 @@ SpBase<elem_type,derived>::is_finite() const
    
    while(it != it_end)
      {
-      if(arma_isfinite(*it) == false)  { return false; }
+      if(arma_isnonfinite(*it))  { return false; }
      ++it;
      }
    }
@@ -850,7 +850,7 @@ SpBase<elem_type,derived>::has_nonfinite() const
    
    while(it != it_end)
      {
-      if(arma_isfinite(*it) == false)  { return true; }
+      if(arma_isnonfinite(*it))  { return true; }
      ++it;
      }
    }
@@ -612,33 +612,26 @@ SpMat<eT>::operator*=(const eT val)
  {
  arma_debug_sigprint();
  
-  if(val != eT(0))
+  sync_csc();
+  invalidate_cache();
+  
+  const uword n_nz = n_nonzero;
+  
+  eT* vals = access::rwp(values);
+  
+  bool has_zero = false;
+  
+  for(uword i=0; i<n_nz; ++i)
    {
-    sync_csc();
-    invalidate_cache();
+    eT& vals_i = vals[i];
    
-    const uword n_nz = n_nonzero;
+    vals_i *= val;
    
-    eT* vals = access::rwp(values);
-    
-    bool has_zero = false;
-    
-    for(uword i=0; i<n_nz; ++i)
-      {
-      eT& vals_i = vals[i];
-      
-      vals_i *= val;
-      
-      if(vals_i == eT(0))  { has_zero = true; }
-      }
-    
-    if(has_zero)  { remove_zeros(); }
-    }
-  else
-    {
-    (*this).zeros();
+    if(vals_i == eT(0))  { has_zero = true; }
    }
  
+  if(has_zero)  { remove_zeros(); }
+  
  return *this;
  }

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


@@ -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)
@@ -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;
@@ -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);
@@ -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();
@@ -36,6 +36,9 @@ class SpToDGlue : public Base< typename T1::elem_type, SpToDGlue<T1, T2, glue_ty
  inline explicit SpToDGlue(const T1& in_A, const T2& in_B);
  inline         ~SpToDGlue();
  
+  template<typename eT2>
+  constexpr bool is_alias(const Mat<eT2>&) const { return false; }
+  
  const T1& A;  //!< first operand;  must be derived from Base or SpBase
  const T2& B;  //!< second operand; must be derived from Base or SpBase
  };
@@ -39,10 +39,13 @@ class SpToDOp : public Base< typename T1::elem_type, SpToDOp<T1, op_type> >
  inline          SpToDOp(const T1& in_m, const uword     in_aux_uword_a, const uword in_aux_uword_b);
  inline         ~SpToDOp();
  
-  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
+  template<typename eT2>
+  constexpr bool is_alias(const Mat<eT2>&) const { return false; }
+  
+  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
  };


@@ -148,24 +148,17 @@ inline
 SpValProxy<T1>&
 SpValProxy<T1>::operator*=(const eT rhs)
  {
-  if(rhs != eT(0))
+  if(val_ptr)
    {
-    if(val_ptr)
-      {
-      // The value already exists and merely needs to be updated.
-      *val_ptr *= rhs;
-      parent.invalidate_cache();
-      check_zero();
-      }
+    *val_ptr *= rhs;
+    parent.invalidate_cache();
+    check_zero();
    }
  else
    {
-    if(val_ptr)
-      {
-      // Since we are multiplying by zero, the value can be deleted.
-      parent.delete_element(row, col);
-      val_ptr = nullptr;
-      }
+    const eT val = eT(0) * rhs;  // in case rhs is inf or nan
+    
+    if(val != eT(0))  { val_ptr = &parent.insert_element(row, col, val); }
    }
  
  return *this;
@@ -178,37 +171,17 @@ inline
 SpValProxy<T1>&
 SpValProxy<T1>::operator/=(const eT rhs)
  {
-  if(rhs != eT(0)) // I hope this is true!
+  if(val_ptr)
    {
-    if(val_ptr)
-      {
-      *val_ptr /= rhs;
-      parent.invalidate_cache();
-      check_zero();
-      }
+    *val_ptr /= rhs;
+    parent.invalidate_cache();
+    check_zero();
    }
  else
    {
-    if(val_ptr)
-      {
-      *val_ptr /= rhs; // That is where it gets ugly.
-      // Now check if it's 0.
-      if(*val_ptr == eT(0))
-        {
-        parent.delete_element(row, col);
-        val_ptr = nullptr;
-        }
-      }
-    else
-      {
-      eT val = eT(0) / rhs; // This may vary depending on type and implementation.
-      
-      if(val != eT(0))
-        {
-        // Ok, now we have to insert it.
-        val_ptr = &parent.insert_element(row, col, val);
-        }
-      }
+    const eT val = eT(0) / rhs;  // in case rhs is zero or nan
+    
+    if(val != eT(0))  { val_ptr = &parent.insert_element(row, col, val); }
    }
  
  return *this;
@@ -65,6 +65,48 @@ arma_isfinite(const std::complex<T>& x)
  }


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

 //
 // wrappers for isinf
@@ -65,6 +65,13 @@ struct arma_config
  #endif
  
  
+  #if defined(ARMA_OPTIMISE_POWEXPR)
+    static constexpr bool optimise_powexpr = true;
+  #else
+    static constexpr bool optimise_powexpr = false;
+  #endif
+  
+  
  #if defined(ARMA_CHECK_CONFORMANCE)
    static constexpr bool check_conform = true;
  #else
@@ -91,6 +91,9 @@ class op_diagmat;
 class op_trimat;
 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;
@@ -128,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;
@@ -358,6 +362,7 @@ struct arma_zeros_indicator   : public arma_initmode_indicator<true > {};
 struct arma_nozeros_indicator : public arma_initmode_indicator<false> {};


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

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


 #define ARMA_VERSION_MAJOR 14
-#define ARMA_VERSION_MINOR 0
+#define ARMA_VERSION_MINOR 6
 #define ARMA_VERSION_PATCH 2
-#define ARMA_VERSION_NAME  "Stochastic Parrot"
+#define ARMA_VERSION_NAME  "Caffe Mocha"



@@ -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;
@@ -43,6 +43,18 @@ class auxlib
  template<typename eT>
  inline static bool inv_tr_rcond(Mat<eT>& A, typename get_pod_type<eT>::result& out_rcond, const uword layout);
  
+  template<typename eT>
+  inline static bool inv_sym(Mat<eT>& A);
+  
+  template<typename T>
+  inline static bool inv_sym(Mat< std::complex<T> >& A);
+  
+  template<typename eT>
+  inline static bool inv_sym_rcond(Mat<eT>& A, eT& out_rcond);
+  
+  template<typename T>
+  inline static bool inv_sym_rcond(Mat< std::complex<T> >& A, T& out_rcond);
+  
  template<typename eT>
  inline static bool inv_sympd(Mat<eT>& A, bool& out_sympd_state);
  
@@ -50,10 +62,10 @@ class auxlib
  inline static bool inv_sympd(Mat<eT>& out, const Mat<eT>& X);
  
  template<typename eT>
-  inline static bool inv_sympd_rcond(Mat<eT>& A, bool& out_sympd_state, eT& out_rcond);
+  inline static bool inv_sympd_rcond(Mat<eT>& A, eT& out_rcond);
  
  template<typename T>
-  inline static bool inv_sympd_rcond(Mat< std::complex<T> >& A, bool& out_sympd_state, T& out_rcond);
+  inline static bool inv_sympd_rcond(Mat< std::complex<T> >& A, T& out_rcond);
  
  
  //
@@ -269,6 +281,20 @@ class auxlib
  
  //
  
+  template<typename T1>
+  inline static bool solve_sym_fast(Mat<typename T1::pod_type>& out, Mat<typename T1::pod_type>& A, const Base<typename T1::pod_type,T1>& B_expr);
+  
+  template<typename T1>
+  inline static bool solve_sym_fast(Mat< std::complex<typename T1::pod_type> >& out, Mat< std::complex<typename T1::pod_type> >& A, const Base< std::complex<typename T1::pod_type>, T1 >& B_expr);
+  
+  template<typename T1>
+  inline static bool solve_sym_rcond(Mat<typename T1::pod_type>& out, typename T1::pod_type& out_rcond, Mat<typename T1::pod_type>& A, const Base<typename T1::pod_type,T1>& B_expr);
+  
+  template<typename T1>
+  inline static bool solve_sym_rcond(Mat< std::complex<typename T1::pod_type> >& out, typename T1::pod_type& out_rcond, Mat< std::complex<typename T1::pod_type> >& A, const Base< std::complex<typename T1::pod_type>,T1>& B_expr);
+  
+  //
+  
  template<typename T1>
  inline static bool solve_sympd_fast(Mat<typename T1::elem_type>& out, Mat<typename T1::elem_type>& A, const Base<typename T1::elem_type,T1>& B_expr);
  
@@ -365,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);
  
  
  //
@@ -378,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
  
@@ -388,10 +421,10 @@ class auxlib
  inline static  T rcond(Mat< std::complex<T> >& A);
  
  template<typename eT>
-  inline static eT rcond_sympd(Mat<eT>& A, bool& calc_ok);
+  inline static eT rcond_sym(Mat<eT>& A);
  
-  template<typename T>
-  inline static  T rcond_sympd(Mat< std::complex<T> >& A, bool& calc_ok);
+  template<typename  T>
+  inline static  T rcond_sym(Mat< std::complex<T> >& A);
  
  template<typename eT>
  inline static eT rcond_trimat(const Mat<eT>& A, const uword layout);
@@ -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
  
@@ -83,10 +83,8 @@
 #endif


-#if (!defined(ARMA_HAVE_CXX14))
-  #if (!defined(ARMA_IGNORE_DEPRECATED_MARKER)) || defined(ARMA_DONT_IGNORE_DEPRECATED_MARKER) || defined(ARMA_DEBUG)
-    #pragma message ("INFO: support for C++11 is deprecated")
-  #endif
+#if (!defined(ARMA_HAVE_CXX14)) && (!defined(ARMA_IGNORE_DEPRECATED_MARKER))
+  #pragma message ("INFO: support for C++11 is deprecated; minimum recommended standard is C++14")
 #endif


@@ -168,6 +168,10 @@
  // gcc 6.1 has proper C++14 support and fixes an OpenMP related bug:
  // 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")
+  #endif
+  
  #define ARMA_GOOD_COMPILER
  
  #undef  arma_hot
@@ -310,7 +314,7 @@

 #if defined(_MSC_VER)
  
-  #if (_MSC_VER < 1900)
+  #if (_MSC_VER < 1910)
    #error "*** newer compiler required ***"
  #endif
  
@@ -473,11 +477,7 @@
 #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) && (!defined(ARMA_DONT_IGNORE_DEPRECATED_MARKER)) && (!defined(ARMA_DEBUG))
+#if defined(ARMA_IGNORE_DEPRECATED_MARKER)
  #undef  arma_deprecated
  #define arma_deprecated

@@ -102,11 +102,8 @@
 //// Uncomment the above line if your BLAS and LAPACK libraries have function names with a trailing underscore.
 //// Conversely, comment it out if the function names don't have a trailing underscore.

-// #define ARMA_BLAS_LONG
-//// Uncomment the above line if your BLAS and LAPACK libraries use "long" instead of "int"
-
 // #define ARMA_BLAS_LONG_LONG
-//// Uncomment the above line if your BLAS and LAPACK libraries use "long long" instead of "int"
+//// Uncomment the above line if your BLAS and LAPACK libraries use 64 bit integers, ie. "long long" instead of "int"

 // #define ARMA_BLAS_NOEXCEPT
 //// Uncomment the above line if you require BLAS functions to have the 'noexcept' specification
@@ -131,7 +128,7 @@
 //// Uncomment the above line to use Intel MKL types for complex numbers.
 //// You will need to include appropriate MKL headers before the Armadillo header.
 //// You may also need to enable or disable the following options:
-//// ARMA_BLAS_LONG, ARMA_BLAS_LONG_LONG, ARMA_USE_FORTRAN_HIDDEN_ARGS
+//// ARMA_BLAS_LONG_LONG, ARMA_USE_FORTRAN_HIDDEN_ARGS

 #if !defined(ARMA_USE_OPENMP)
 // #define ARMA_USE_OPENMP
@@ -169,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.
@@ -214,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

@@ -293,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
@@ -306,6 +308,24 @@
  #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
+  #define ARMA_BLAS_64BIT_INT
+#endif
+
 #if defined(ARMA_DONT_OPTIMISE_BAND) || defined(ARMA_DONT_OPTIMISE_SOLVE_BAND)
  #undef ARMA_OPTIMISE_BAND
 #endif
@@ -318,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")
@@ -330,6 +354,10 @@
  #undef ARMA_CHECK_NONFINITE
 #endif

+#if defined(ARMA_DONT_IGNORE_DEPRECATED_MARKER)
+  #undef ARMA_IGNORE_DEPRECATED_MARKER
+#endif
+
 #if defined(ARMA_NO_DEBUG)
  #undef ARMA_DEBUG
  #undef ARMA_EXTRA_DEBUG
@@ -347,6 +375,8 @@
  
  #undef  ARMA_WARN_LEVEL
  #define ARMA_WARN_LEVEL 3
+  
+  #undef  ARMA_IGNORE_DEPRECATED_MARKER
 #endif

 #if defined(ARMA_DONT_PRINT_EXCEPTIONS)
@@ -357,6 +387,15 @@
  #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 ("WARNING: option ARMA_CRIPPLED_LAPACK is deprecated and will be removed")
+  #endif
+#endif
+

 // if Armadillo was installed on this system via CMake and ARMA_USE_WRAPPER is not defined,
 // ARMA_AUX_LIBS lists the libraries required by Armadillo on this system, and
@@ -102,11 +102,8 @@
 //// Uncomment the above line if your BLAS and LAPACK libraries have function names with a trailing underscore.
 //// Conversely, comment it out if the function names don't have a trailing underscore.

-// #define ARMA_BLAS_LONG
-//// Uncomment the above line if your BLAS and LAPACK libraries use "long" instead of "int"
-
 // #define ARMA_BLAS_LONG_LONG
-//// Uncomment the above line if your BLAS and LAPACK libraries use "long long" instead of "int"
+//// Uncomment the above line if your BLAS and LAPACK libraries use 64 bit integers, ie. "long long" instead of "int"

 // #define ARMA_BLAS_NOEXCEPT
 //// Uncomment the above line if you require BLAS functions to have the 'noexcept' specification
@@ -131,7 +128,7 @@
 //// Uncomment the above line to use Intel MKL types for complex numbers.
 //// You will need to include appropriate MKL headers before the Armadillo header.
 //// You may also need to enable or disable the following options:
-//// ARMA_BLAS_LONG, ARMA_BLAS_LONG_LONG, ARMA_USE_FORTRAN_HIDDEN_ARGS
+//// ARMA_BLAS_LONG_LONG, ARMA_USE_FORTRAN_HIDDEN_ARGS

 #if !defined(ARMA_USE_OPENMP)
 // #define ARMA_USE_OPENMP
@@ -169,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.
@@ -214,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

@@ -293,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
@@ -306,6 +308,24 @@
  #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
+  #define ARMA_BLAS_64BIT_INT
+#endif
+
 #if defined(ARMA_DONT_OPTIMISE_BAND) || defined(ARMA_DONT_OPTIMISE_SOLVE_BAND)
  #undef ARMA_OPTIMISE_BAND
 #endif
@@ -318,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")
@@ -330,6 +354,10 @@
  #undef ARMA_CHECK_NONFINITE
 #endif

+#if defined(ARMA_DONT_IGNORE_DEPRECATED_MARKER)
+  #undef ARMA_IGNORE_DEPRECATED_MARKER
+#endif
+
 #if defined(ARMA_NO_DEBUG)
  #undef ARMA_DEBUG
  #undef ARMA_EXTRA_DEBUG
@@ -347,6 +375,8 @@
  
  #undef  ARMA_WARN_LEVEL
  #define ARMA_WARN_LEVEL 3
+  
+  #undef  ARMA_IGNORE_DEPRECATED_MARKER
 #endif

 #if defined(ARMA_DONT_PRINT_EXCEPTIONS)
@@ -357,6 +387,15 @@
  #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 ("WARNING: option ARMA_CRIPPLED_LAPACK is deprecated and will be removed")
+  #endif
+#endif
+

 // if Armadillo was installed on this system via CMake and ARMA_USE_WRAPPER is not defined,
 // ARMA_AUX_LIBS lists the libraries required by Armadillo on this system, and
@@ -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;
@@ -269,6 +269,35 @@
  #define arma_cpstrf cpstrf
  #define arma_zpstrf zpstrf
  
+  #define arma_ssytrf ssytrf
+  #define arma_dsytrf dsytrf
+  
+  #define arma_chetrf chetrf
+  #define arma_zhetrf zhetrf
+  
+  #define arma_ssytrs ssytrs
+  #define arma_dsytrs dsytrs
+  
+  #define arma_chetrs chetrs
+  #define arma_zhetrs zhetrs
+  
+  #define arma_ssytri ssytri
+  #define arma_dsytri dsytri
+  
+  #define arma_chetri chetri
+  #define arma_zhetri zhetri
+  
+  #define arma_ssycon ssycon
+  #define arma_dsycon dsycon
+  
+  #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
@@ -505,6 +534,35 @@
  #define arma_cpstrf CPSTRF
  #define arma_zpstrf ZPSTRF
  
+  #define arma_ssytrf SSYTRF
+  #define arma_dsytrf DSYTRF
+  
+  #define arma_chetrf CHETRF
+  #define arma_zhetrf ZHETRF
+  
+  #define arma_ssytrs SSYTRS
+  #define arma_dsytrs DSYTRS
+  
+  #define arma_chetrs CHETRS
+  #define arma_zhetrs ZHETRS
+  
+  #define arma_ssytri SSYTRI
+  #define arma_dsytri DSYTRI
+  
+  #define arma_chetri CHETRI
+  #define arma_zhetri ZHETRI
+  
+  #define arma_ssycon SSYCON
+  #define arma_dsycon DSYCON
+  
+  #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


@@ -846,6 +904,44 @@ extern "C"
  void arma_fortran(arma_cpstrf)(const char* uplo, const blas_int* n, blas_cxf* a, const blas_int* lda, blas_int* piv, blas_int* rank, const  float* tol,  float* work, blas_int* info, blas_len uplo_len) ARMA_NOEXCEPT;
  void arma_fortran(arma_zpstrf)(const char* uplo, const blas_int* n, blas_cxd* a, const blas_int* lda, blas_int* piv, blas_int* rank, const double* tol, double* work, blas_int* info, blas_len uplo_len) ARMA_NOEXCEPT;
  
+  // factorisation of symmetric matrix (real)
+  void arma_fortran(arma_ssytrf)(const char* uplo, const blas_int* n, float*    a, const blas_int* lda, blas_int* ipiv, float*    work, const blas_int* lwork, blas_int* info, blas_len uplo_len) ARMA_NOEXCEPT;
+  void arma_fortran(arma_dsytrf)(const char* uplo, const blas_int* n, double*   a, const blas_int* lda, blas_int* ipiv, double*   work, const blas_int* lwork, blas_int* info, blas_len uplo_len) ARMA_NOEXCEPT;
+  
+  // factorisation of hermitian matrix (complex)
+  void arma_fortran(arma_chetrf)(const char* uplo, const blas_int* n, blas_cxf* a, const blas_int* lda, blas_int* ipiv, blas_cxf* work, const blas_int* lwork, blas_int* info, blas_len uplo_len) ARMA_NOEXCEPT;
+  void arma_fortran(arma_zhetrf)(const char* uplo, const blas_int* n, blas_cxd* a, const blas_int* lda, blas_int* ipiv, blas_cxd* work, const blas_int* lwork, blas_int* info, blas_len uplo_len) ARMA_NOEXCEPT;
+  
+  // solve system using pre-computed factorisation (real)
+  void arma_fortran(arma_ssytrs)(const char* uplo, const blas_int* n, const blas_int* nrhs, const float*    a, const blas_int* lda, const blas_int* ipiv, float*  b, const blas_int* ldb, blas_int* info, blas_len uplo_len) ARMA_NOEXCEPT;
+  void arma_fortran(arma_dsytrs)(const char* uplo, const blas_int* n, const blas_int* nrhs, const double*   a, const blas_int* lda, const blas_int* ipiv, double* b, const blas_int* ldb, blas_int* info, blas_len uplo_len) ARMA_NOEXCEPT;
+  
+  // solve system using pre-computed factorisation (complex)
+  void arma_fortran(arma_chetrs)(const char* uplo, const blas_int* n, const blas_int* nrhs, const blas_cxf* a, const blas_int* lda, const blas_int* ipiv, blas_cxf* b, const blas_int* ldb, blas_int* info, blas_len uplo_len) ARMA_NOEXCEPT;
+  void arma_fortran(arma_zhetrs)(const char* uplo, const blas_int* n, const blas_int* nrhs, const blas_cxd* a, const blas_int* lda, const blas_int* ipiv, blas_cxd* b, const blas_int* ldb, blas_int* info, blas_len uplo_len) ARMA_NOEXCEPT;
+  
+  // inverse of symmetric matrix using pre-computed factorisation (real)
+  void arma_fortran(arma_ssytri)(const char* uplo, const blas_int* n, float*    a, const blas_int* lda, const blas_int* ipiv, float*    work, blas_int* info, blas_len uplo_len) ARMA_NOEXCEPT;
+  void arma_fortran(arma_dsytri)(const char* uplo, const blas_int* n, double*   a, const blas_int* lda, const blas_int* ipiv, double*   work, blas_int* info, blas_len uplo_len) ARMA_NOEXCEPT;
+  
+  // inverse of hermitian matrix using pre-computed factorisation (complex)
+  void arma_fortran(arma_chetri)(const char* uplo, const blas_int* n, blas_cxf* a, const blas_int* lda, const blas_int* ipiv, blas_cxf* work, blas_int* info, blas_len uplo_len) ARMA_NOEXCEPT;
+  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, 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, const blas_cxf* a, const blas_int* lda, const blas_int* ipiv, const float*  anorm, float*  rcond, blas_cxf* work, blas_int* info, blas_len uplo_len) ARMA_NOEXCEPT;
+  void arma_fortran(arma_zhecon)(const char* uplo, const blas_int* n, const blas_cxd* a, const blas_int* lda, const blas_int* ipiv, const double* anorm, double* rcond, blas_cxd* work, blas_int* info, blas_len uplo_len) ARMA_NOEXCEPT;
+  
+  // matrix balance
+  void arma_fortran(arma_sgebal)(const char* job, const blas_int* n,    float* a, const blas_int* lda, blas_int* ilo, blas_int* ihi,  float* scale, blas_int* info, blas_len job_len) ARMA_NOEXCEPT;
+  void arma_fortran(arma_dgebal)(const char* job, const blas_int* n,   double* a, const blas_int* lda, blas_int* ilo, blas_int* ihi, double* scale, blas_int* info, blas_len job_len) ARMA_NOEXCEPT;
+  void arma_fortran(arma_cgebal)(const char* job, const blas_int* n, blas_cxf* a, const blas_int* lda, blas_int* ilo, blas_int* ihi,  float* scale, blas_int* info, blas_len job_len) ARMA_NOEXCEPT;
+  void arma_fortran(arma_zgebal)(const char* job, const blas_int* n, blas_cxd* a, const blas_int* lda, blas_int* ilo, blas_int* ihi, double* scale, blas_int* info, blas_len job_len) ARMA_NOEXCEPT;
+  
 #else
  
  // prototypes without hidden arguments
@@ -1170,6 +1266,44 @@ extern "C"
  void arma_fortran(arma_cpstrf)(const char* uplo, const blas_int* n, blas_cxf* a, const blas_int* lda, blas_int* piv, blas_int* rank, const  float* tol,  float* work, blas_int* info) ARMA_NOEXCEPT;
  void arma_fortran(arma_zpstrf)(const char* uplo, const blas_int* n, blas_cxd* a, const blas_int* lda, blas_int* piv, blas_int* rank, const double* tol, double* work, blas_int* info) ARMA_NOEXCEPT;
  
+  // factorisation of symmetric matrix (real)
+  void arma_fortran(arma_ssytrf)(const char* uplo, const blas_int* n, float*    a, const blas_int* lda, blas_int* ipiv, float*    work, const blas_int* lwork, blas_int* info) ARMA_NOEXCEPT;
+  void arma_fortran(arma_dsytrf)(const char* uplo, const blas_int* n, double*   a, const blas_int* lda, blas_int* ipiv, double*   work, const blas_int* lwork, blas_int* info) ARMA_NOEXCEPT;
+  
+  // factorisation of hermitian matrix (complex)
+  void arma_fortran(arma_chetrf)(const char* uplo, const blas_int* n, blas_cxf* a, const blas_int* lda, blas_int* ipiv, blas_cxf* work, const blas_int* lwork, blas_int* info) ARMA_NOEXCEPT;
+  void arma_fortran(arma_zhetrf)(const char* uplo, const blas_int* n, blas_cxd* a, const blas_int* lda, blas_int* ipiv, blas_cxd* work, const blas_int* lwork, blas_int* info) ARMA_NOEXCEPT;
+  
+  // solve system using pre-computed factorisation (real)
+  void arma_fortran(arma_ssytrs)(const char* uplo, const blas_int* n, const blas_int* nrhs, const float*    a, const blas_int* lda, const blas_int* ipiv, float*  b, const blas_int* ldb, blas_int* info) ARMA_NOEXCEPT;
+  void arma_fortran(arma_dsytrs)(const char* uplo, const blas_int* n, const blas_int* nrhs, const double*   a, const blas_int* lda, const blas_int* ipiv, double* b, const blas_int* ldb, blas_int* info) ARMA_NOEXCEPT;
+  
+  // solve system using pre-computed factorisation (complex)
+  void arma_fortran(arma_zhetrs)(const char* uplo, const blas_int* n, const blas_int* nrhs, const blas_cxd* a, const blas_int* lda, const blas_int* ipiv, blas_cxd* b, const blas_int* ldb, blas_int* info) ARMA_NOEXCEPT;
+  void arma_fortran(arma_chetrs)(const char* uplo, const blas_int* n, const blas_int* nrhs, const blas_cxf* a, const blas_int* lda, const blas_int* ipiv, blas_cxf* b, const blas_int* ldb, blas_int* info) ARMA_NOEXCEPT;
+  
+  // inverse of symmetric matrix using pre-computed factorisation (real)
+  void arma_fortran(arma_ssytri)(const char* uplo, const blas_int* n, float*    a, const blas_int* lda, const blas_int* ipiv, float*    work, blas_int* info) ARMA_NOEXCEPT;
+  void arma_fortran(arma_dsytri)(const char* uplo, const blas_int* n, double*   a, const blas_int* lda, const blas_int* ipiv, double*   work, blas_int* info) ARMA_NOEXCEPT;
+  
+  // inverse of hermitian matrix using pre-computed factorisation (complex)
+  void arma_fortran(arma_chetri)(const char* uplo, const blas_int* n, blas_cxf* a, const blas_int* lda, const blas_int* ipiv, blas_cxf* work, blas_int* info) ARMA_NOEXCEPT;
+  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, 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, const blas_cxf* a, const blas_int* lda, const blas_int* ipiv, const float*  anorm, float*  rcond, blas_cxf* work, blas_int* info) ARMA_NOEXCEPT;
+  void arma_fortran(arma_zhecon)(const char* uplo, const blas_int* n, const blas_cxd* a, const blas_int* lda, const blas_int* ipiv, const double* anorm, double* rcond, blas_cxd* work, blas_int* info) ARMA_NOEXCEPT;
+  
+  // matrix balance
+  void arma_fortran(arma_sgebal)(const char* job, const blas_int* n,    float* a, const blas_int* lda, blas_int* ilo, blas_int* ihi,  float* scale, blas_int* info) ARMA_NOEXCEPT;
+  void arma_fortran(arma_dgebal)(const char* job, const blas_int* n,   double* a, const blas_int* lda, blas_int* ilo, blas_int* ihi, double* scale, blas_int* info) ARMA_NOEXCEPT;
+  void arma_fortran(arma_cgebal)(const char* job, const blas_int* n, blas_cxf* a, const blas_int* lda, blas_int* ilo, blas_int* ihi,  float* scale, blas_int* info) ARMA_NOEXCEPT;
+  void arma_fortran(arma_zgebal)(const char* job, const blas_int* n, blas_cxd* a, const blas_int* lda, blas_int* ilo, blas_int* ihi, double* scale, blas_int* info) ARMA_NOEXCEPT;
+  
 #endif
 }

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

 extern "C"
  {
-  extern void arma_wrapper(sgssv)(superlu::superlu_options_t*, superlu::SuperMatrix*, int*, int*, superlu::SuperMatrix*, superlu::SuperMatrix*, superlu::SuperMatrix*, superlu::SuperLUStat_t*, int*);
-  extern void arma_wrapper(dgssv)(superlu::superlu_options_t*, superlu::SuperMatrix*, int*, int*, superlu::SuperMatrix*, superlu::SuperMatrix*, superlu::SuperMatrix*, superlu::SuperLUStat_t*, int*);
-  extern void arma_wrapper(cgssv)(superlu::superlu_options_t*, superlu::SuperMatrix*, int*, int*, superlu::SuperMatrix*, superlu::SuperMatrix*, superlu::SuperMatrix*, superlu::SuperLUStat_t*, int*);
-  extern void arma_wrapper(zgssv)(superlu::superlu_options_t*, superlu::SuperMatrix*, int*, int*, superlu::SuperMatrix*, superlu::SuperMatrix*, superlu::SuperMatrix*, superlu::SuperLUStat_t*, int*);
+  extern void arma_wrapper(sgssv)(superlu::superlu_options_t*, superlu::SuperMatrix*, int*, int*, superlu::SuperMatrix*, superlu::SuperMatrix*, superlu::SuperMatrix*, superlu::SuperLUStat_t*, superlu::int_t*);
+  extern void arma_wrapper(dgssv)(superlu::superlu_options_t*, superlu::SuperMatrix*, int*, int*, superlu::SuperMatrix*, superlu::SuperMatrix*, superlu::SuperMatrix*, superlu::SuperLUStat_t*, superlu::int_t*);
+  extern void arma_wrapper(cgssv)(superlu::superlu_options_t*, superlu::SuperMatrix*, int*, int*, superlu::SuperMatrix*, superlu::SuperMatrix*, superlu::SuperMatrix*, superlu::SuperLUStat_t*, superlu::int_t*);
+  extern void arma_wrapper(zgssv)(superlu::superlu_options_t*, superlu::SuperMatrix*, int*, int*, superlu::SuperMatrix*, superlu::SuperMatrix*, superlu::SuperMatrix*, superlu::SuperLUStat_t*, superlu::int_t*);
  
-  extern void arma_wrapper(sgssvx)(superlu::superlu_options_t*, superlu::SuperMatrix*, int*, int*, int*, char*,  float*,  float*, superlu::SuperMatrix*, superlu::SuperMatrix*, void*, int, superlu::SuperMatrix*, superlu::SuperMatrix*,  float*,  float*,  float*,  float*, superlu::GlobalLU_t*, superlu::mem_usage_t*, superlu::SuperLUStat_t*, int*);
-  extern void arma_wrapper(dgssvx)(superlu::superlu_options_t*, superlu::SuperMatrix*, int*, int*, int*, char*, double*, double*, superlu::SuperMatrix*, superlu::SuperMatrix*, void*, int, superlu::SuperMatrix*, superlu::SuperMatrix*, double*, double*, double*, double*, superlu::GlobalLU_t*, superlu::mem_usage_t*, superlu::SuperLUStat_t*, int*);
-  extern void arma_wrapper(cgssvx)(superlu::superlu_options_t*, superlu::SuperMatrix*, int*, int*, int*, char*,  float*,  float*, superlu::SuperMatrix*, superlu::SuperMatrix*, void*, int, superlu::SuperMatrix*, superlu::SuperMatrix*,  float*,  float*,  float*,  float*, superlu::GlobalLU_t*, superlu::mem_usage_t*, superlu::SuperLUStat_t*, int*);
-  extern void arma_wrapper(zgssvx)(superlu::superlu_options_t*, superlu::SuperMatrix*, int*, int*, int*, char*, double*, double*, superlu::SuperMatrix*, superlu::SuperMatrix*, void*, int, superlu::SuperMatrix*, superlu::SuperMatrix*, double*, double*, double*, double*, superlu::GlobalLU_t*, superlu::mem_usage_t*, superlu::SuperLUStat_t*, int*);
+  extern void arma_wrapper(sgssvx)(superlu::superlu_options_t*, superlu::SuperMatrix*, int*, int*, int*, char*,  float*,  float*, superlu::SuperMatrix*, superlu::SuperMatrix*, void*, superlu::int_t, superlu::SuperMatrix*, superlu::SuperMatrix*,  float*,  float*,  float*,  float*, superlu::GlobalLU_t*, superlu::mem_usage_t*, superlu::SuperLUStat_t*, superlu::int_t*);
+  extern void arma_wrapper(dgssvx)(superlu::superlu_options_t*, superlu::SuperMatrix*, int*, int*, int*, char*, double*, double*, superlu::SuperMatrix*, superlu::SuperMatrix*, void*, superlu::int_t, superlu::SuperMatrix*, superlu::SuperMatrix*, double*, double*, double*, double*, superlu::GlobalLU_t*, superlu::mem_usage_t*, superlu::SuperLUStat_t*, superlu::int_t*);
+  extern void arma_wrapper(cgssvx)(superlu::superlu_options_t*, superlu::SuperMatrix*, int*, int*, int*, char*,  float*,  float*, superlu::SuperMatrix*, superlu::SuperMatrix*, void*, superlu::int_t, superlu::SuperMatrix*, superlu::SuperMatrix*,  float*,  float*,  float*,  float*, superlu::GlobalLU_t*, superlu::mem_usage_t*, superlu::SuperLUStat_t*, superlu::int_t*);
+  extern void arma_wrapper(zgssvx)(superlu::superlu_options_t*, superlu::SuperMatrix*, int*, int*, int*, char*, double*, double*, superlu::SuperMatrix*, superlu::SuperMatrix*, void*, superlu::int_t, superlu::SuperMatrix*, superlu::SuperMatrix*, double*, double*, double*, double*, superlu::GlobalLU_t*, superlu::mem_usage_t*, superlu::SuperLUStat_t*, superlu::int_t*);
  
-  extern void arma_wrapper(sgstrf)(superlu::superlu_options_t*, superlu::SuperMatrix*, int, int, int*, void*, int, int*, int*, superlu::SuperMatrix*, superlu::SuperMatrix*, superlu::GlobalLU_t*, superlu::SuperLUStat_t*, int*);
-  extern void arma_wrapper(dgstrf)(superlu::superlu_options_t*, superlu::SuperMatrix*, int, int, int*, void*, int, int*, int*, superlu::SuperMatrix*, superlu::SuperMatrix*, superlu::GlobalLU_t*, superlu::SuperLUStat_t*, int*);
-  extern void arma_wrapper(cgstrf)(superlu::superlu_options_t*, superlu::SuperMatrix*, int, int, int*, void*, int, int*, int*, superlu::SuperMatrix*, superlu::SuperMatrix*, superlu::GlobalLU_t*, superlu::SuperLUStat_t*, int*);
-  extern void arma_wrapper(zgstrf)(superlu::superlu_options_t*, superlu::SuperMatrix*, int, int, int*, void*, int, int*, int*, superlu::SuperMatrix*, superlu::SuperMatrix*, superlu::GlobalLU_t*, superlu::SuperLUStat_t*, int*);
+  extern void arma_wrapper(sgstrf)(superlu::superlu_options_t*, superlu::SuperMatrix*, int, int, int*, void*, superlu::int_t, int*, int*, superlu::SuperMatrix*, superlu::SuperMatrix*, superlu::GlobalLU_t*, superlu::SuperLUStat_t*, superlu::int_t*);
+  extern void arma_wrapper(dgstrf)(superlu::superlu_options_t*, superlu::SuperMatrix*, int, int, int*, void*, superlu::int_t, int*, int*, superlu::SuperMatrix*, superlu::SuperMatrix*, superlu::GlobalLU_t*, superlu::SuperLUStat_t*, superlu::int_t*);
+  extern void arma_wrapper(cgstrf)(superlu::superlu_options_t*, superlu::SuperMatrix*, int, int, int*, void*, superlu::int_t, int*, int*, superlu::SuperMatrix*, superlu::SuperMatrix*, superlu::GlobalLU_t*, superlu::SuperLUStat_t*, superlu::int_t*);
+  extern void arma_wrapper(zgstrf)(superlu::superlu_options_t*, superlu::SuperMatrix*, int, int, int*, void*, superlu::int_t, int*, int*, superlu::SuperMatrix*, superlu::SuperMatrix*, superlu::GlobalLU_t*, superlu::SuperLUStat_t*, superlu::int_t*);
  
  extern void arma_wrapper(sgstrs)(superlu::trans_t, superlu::SuperMatrix*, superlu::SuperMatrix*, int*, int*, superlu::SuperMatrix*, superlu::SuperLUStat_t*, int*);
  extern void arma_wrapper(dgstrs)(superlu::trans_t, superlu::SuperMatrix*, superlu::SuperMatrix*, int*, int*, superlu::SuperMatrix*, superlu::SuperLUStat_t*, int*);
@@ -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;
@@ -108,6 +108,9 @@ class diagview : public Base< eT, diagview<eT> >
  inline static void schur_inplace(Mat<eT>& out, const diagview& in);
  inline static void   div_inplace(Mat<eT>& out, const diagview& in);
  
+  template<typename eT2>
+  inline bool is_alias(const Mat<eT2>& X) const;
+  
  
  friend class Mat<eT>;
  friend class subview<eT>;
@@ -236,13 +236,13 @@ diagview<eT>::operator= (const Base<eT,T1>& o)
    "diagview: given object has incompatible size"
    );
  
-  const bool is_alias = P.is_alias(d_m);
+  const bool have_alias = P.is_alias(d_m);
  
-  if(is_alias)  { arma_debug_print("aliasing detected"); }
+  if(have_alias)  { arma_debug_print("aliasing detected"); }
  
-  if( (is_Mat<typename Proxy<T1>::stored_type>::value) || (Proxy<T1>::use_at) || (is_alias) )
+  if( (is_Mat<typename Proxy<T1>::stored_type>::value) || (Proxy<T1>::use_at) || (have_alias) )
    {
-    const unwrap_check<typename Proxy<T1>::stored_type> tmp(P.Q, is_alias);
+    const unwrap_check<typename Proxy<T1>::stored_type> tmp(P.Q, have_alias);
    const Mat<eT>& x = tmp.M;
    
    const eT* x_mem = x.memptr();
@@ -309,13 +309,13 @@ diagview<eT>::operator+=(const Base<eT,T1>& o)
    "diagview: given object has incompatible size"
    );
  
-  const bool is_alias = P.is_alias(d_m);
+  const bool have_alias = P.is_alias(d_m);
  
-  if(is_alias)  { arma_debug_print("aliasing detected"); }
+  if(have_alias)  { arma_debug_print("aliasing detected"); }
  
-  if( (is_Mat<typename Proxy<T1>::stored_type>::value) || (Proxy<T1>::use_at) || (is_alias) )
+  if( (is_Mat<typename Proxy<T1>::stored_type>::value) || (Proxy<T1>::use_at) || (have_alias) )
    {
-    const unwrap_check<typename Proxy<T1>::stored_type> tmp(P.Q, is_alias);
+    const unwrap_check<typename Proxy<T1>::stored_type> tmp(P.Q, have_alias);
    const Mat<eT>& x = tmp.M;
    
    const eT* x_mem = x.memptr();
@@ -382,13 +382,13 @@ diagview<eT>::operator-=(const Base<eT,T1>& o)
    "diagview: given object has incompatible size"
    );
  
-  const bool is_alias = P.is_alias(d_m);
+  const bool have_alias = P.is_alias(d_m);
  
-  if(is_alias)  { arma_debug_print("aliasing detected"); }
+  if(have_alias)  { arma_debug_print("aliasing detected"); }
  
-  if( (is_Mat<typename Proxy<T1>::stored_type>::value) || (Proxy<T1>::use_at) || (is_alias) )
+  if( (is_Mat<typename Proxy<T1>::stored_type>::value) || (Proxy<T1>::use_at) || (have_alias) )
    {
-    const unwrap_check<typename Proxy<T1>::stored_type> tmp(P.Q, is_alias);
+    const unwrap_check<typename Proxy<T1>::stored_type> tmp(P.Q, have_alias);
    const Mat<eT>& x = tmp.M;
    
    const eT* x_mem = x.memptr();
@@ -455,13 +455,13 @@ diagview<eT>::operator%=(const Base<eT,T1>& o)
    "diagview: given object has incompatible size"
    );
  
-  const bool is_alias = P.is_alias(d_m);
+  const bool have_alias = P.is_alias(d_m);
  
-  if(is_alias)  { arma_debug_print("aliasing detected"); }
+  if(have_alias)  { arma_debug_print("aliasing detected"); }
  
-  if( (is_Mat<typename Proxy<T1>::stored_type>::value) || (Proxy<T1>::use_at) || (is_alias) )
+  if( (is_Mat<typename Proxy<T1>::stored_type>::value) || (Proxy<T1>::use_at) || (have_alias) )
    {
-    const unwrap_check<typename Proxy<T1>::stored_type> tmp(P.Q, is_alias);
+    const unwrap_check<typename Proxy<T1>::stored_type> tmp(P.Q, have_alias);
    const Mat<eT>& x = tmp.M;
    
    const eT* x_mem = x.memptr();
@@ -528,13 +528,13 @@ diagview<eT>::operator/=(const Base<eT,T1>& o)
    "diagview: given object has incompatible size"
    );
  
-  const bool is_alias = P.is_alias(d_m);
+  const bool have_alias = P.is_alias(d_m);
  
-  if(is_alias)  { arma_debug_print("aliasing detected"); }
+  if(have_alias)  { arma_debug_print("aliasing detected"); }
  
-  if( (is_Mat<typename Proxy<T1>::stored_type>::value) || (Proxy<T1>::use_at) || (is_alias) )
+  if( (is_Mat<typename Proxy<T1>::stored_type>::value) || (Proxy<T1>::use_at) || (have_alias) )
    {
-    const unwrap_check<typename Proxy<T1>::stored_type> tmp(P.Q, is_alias);
+    const unwrap_check<typename Proxy<T1>::stored_type> tmp(P.Q, have_alias);
    const Mat<eT>& x = tmp.M;
    
    const eT* x_mem = x.memptr();
@@ -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,12 +1018,29 @@ 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];
    }
  }



+template<typename eT>
+template<typename eT2>
+inline
+bool
+diagview<eT>::is_alias(const Mat<eT2>& X) const
+  {
+  arma_debug_sigprint();
+  
+  return m.is_alias(X);
+  }
+
+
+
 //! @}
@@ -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; }
@@ -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);
@@ -38,12 +38,15 @@ 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);
  
+  template<typename eT2>
+  inline bool is_alias(const Mat<eT2>& X) const;
+  
  arma_inline uword get_n_rows() const;
  arma_inline uword get_n_cols() const;
  arma_inline uword get_n_elem() const;
@@ -49,6 +49,17 @@ eGlue<T1,T2,eglue_type>::eGlue(const T1& in_A, const T2& in_B)



+template<typename T1, typename T2, typename eglue_type>
+template<typename eT2>
+inline
+bool
+eGlue<T1,T2,eglue_type>::is_alias(const Mat<eT2>& X) const
+  {
+  return (P1.is_alias(X) || P2.is_alias(X));
+  }
+
+
+
 template<typename T1, typename T2, typename eglue_type>
 arma_inline
 uword
@@ -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);
@@ -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);
@@ -50,6 +49,9 @@ class eOp : public Base< typename T1::elem_type, eOp<T1, eop_type> >
  inline          eOp(const T1& in_m, const uword in_aux_uword_a, const uword in_aux_uword_b);
  inline          eOp(const T1& in_m, const elem_type in_aux, const uword in_aux_uword_a, const uword in_aux_uword_b);
  
+  template<typename eT2>
+  inline bool is_alias(const Mat<eT2>& X) const;
+  
  arma_inline uword get_n_rows() const;
  arma_inline uword get_n_cols() const;
  arma_inline uword get_n_elem() const;
@@ -74,7 +74,20 @@ eOp<T1, eop_type>::~eOp()
  arma_debug_sigprint();
  }

+
+
+template<typename T1, typename eop_type>
+template<typename eT2>
+inline
+bool
+eOp<T1, eop_type>::is_alias(const Mat<eT2>& X) const
+  {
+  arma_debug_sigprint();
  
+  return P.is_alias(X);
+  }
+
+

 template<typename T1, typename eop_type>
 arma_inline
@@ -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();
@@ -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
    
@@ -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();
@@ -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
@@ -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);
  
@@ -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; }
    }
  }

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


@@ -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;
+  }
+
+
+
 //! @}
@@ -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
      }
    }
@@ -147,50 +147,10 @@ logspace(const double A, const double B, const uword N = 50u)



-//
-// log_exp_add
-
-template<typename eT>
-arma_warn_unused
-inline
-typename arma_real_only<eT>::result
-log_add_exp(eT log_a, eT log_b)
-  {
-  if(log_a < log_b)
-    {
-    std::swap(log_a, log_b);
-    }
-  
-  const eT negdelta = log_b - log_a;
-  
-  if( (negdelta < Datum<eT>::log_min) || (arma_isfinite(negdelta) == false) )
-    {
-    return log_a;
-    }
-  else
-    {
-    return (log_a + std::log1p(std::exp(negdelta)));
-    }
-  }
-
-
-
-// for compatibility with earlier versions
-template<typename eT>
-arma_warn_unused
-inline
-typename arma_real_only<eT>::result
-log_add(eT log_a, eT log_b)
-  {
-  return log_add_exp(log_a, log_b);
-  }
-  
-
-
 //! 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)
  {
@@ -203,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)
@@ -217,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)
@@ -231,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)
@@ -584,4 +544,54 @@ affmul(const T1& A, const T2& B)



+namespace priv
+  {
+  // internal use only
+  template<typename eT>
+  arma_warn_unused
+  inline
+  typename arma_real_only<eT>::result
+  internal_log_add_exp(eT log_a, eT log_b)
+    {
+    if(log_a < log_b)  { std::swap(log_a, log_b); }
+    
+    const eT negdelta = log_b - log_a;
+    
+    if( (negdelta < Datum<eT>::log_min) || arma_isnonfinite(negdelta) )
+      {
+      return log_a;
+      }
+    else
+      {
+      return (log_a + std::log1p(std::exp(negdelta)));
+      }
+    }
+  }
+
+
+
+// DO NOT USE; kept only for compatibility with old user code
+template<typename eT>
+arma_deprecated
+inline
+typename arma_real_only<eT>::result
+log_add_exp(eT log_a, eT log_b)
+  {
+  return priv::internal_log_add_exp(log_a, log_b);
+  }
+
+
+
+// DO NOT USE; kept only for compatibility with old user code
+template<typename eT>
+arma_deprecated
+inline
+typename arma_real_only<eT>::result
+log_add(eT log_a, eT log_b)
+  {
+  return priv::internal_log_add_exp(log_a, log_b);
+  }
+
+
+
 //! @}
@@ -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);
+  }
+
+
+
+//! @}
@@ -23,7 +23,7 @@

 template<typename T1, typename T2>
 arma_warn_unused
-arma_inline
+inline
 typename
 enable_if2
  <
@@ -44,7 +44,7 @@ pow


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

 template<typename T1, typename T2>
 arma_warn_unused
-arma_inline
+inline
 const GlueCube<T1, T2, glue_powext>
 pow
  (
@@ -78,7 +78,7 @@ pow


 template<typename eT, typename T2>
-arma_warn_unused
+arma_deprecated
 inline
 Cube<eT>
 pow
@@ -100,7 +100,7 @@ pow

 template<typename T1, typename T2>
 arma_warn_unused
-arma_inline
+inline
 typename
 enable_if2
  <
@@ -121,7 +121,7 @@ pow


 template<typename parent, unsigned int mode, typename T2>
-arma_warn_unused
+arma_frown("refactor your code to use pow() in conjunction with repmat()")
 inline
 typename
 enable_if2
@@ -144,7 +144,7 @@ pow

 template<typename T1, typename T2>
 arma_warn_unused
-arma_inline
+inline
 const mtGlueCube<typename T1::elem_type, T1, T2, glue_powext_cx>
 pow
  (
@@ -160,7 +160,7 @@ pow


 template<typename T, typename T2>
-arma_warn_unused
+arma_deprecated
 inline
 Cube< std::complex<T> >
 pow
@@ -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;
    }
  
@@ -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);
+  }
+
+
+
 //! @}
@@ -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);
+  }
+
+
+
 //! @}
@@ -21,6 +21,9 @@
 //! @{


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



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



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


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



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


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



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


@@ -139,13 +139,13 @@ glue_mvnrnd::apply_noalias(Mat<eT>& out, const Mat<eT>& M, const Mat<eT>& C, con
    
    const eT tol = eT(-100) * Datum<eT>::eps * norm(C, "fro");
    
-    if(arma_isfinite(tol) == false)  { return false; }
+    if(arma_isnonfinite(tol))  { return false; }
    
    for(uword i=0; i<eigval_n_elem; ++i)
      {
      const eT val = eigval_mem[i];
      
-      if( (val < tol) || (arma_isfinite(val) == false) )  { return false; }
+      if( (val < tol) || arma_isnonfinite(val) )  { return false; }
      }
    
    for(uword i=0; i<eigval_n_elem; ++i)  { if(eigval_mem[i] < eT(0))  { eigval_mem[i] = eT(0); } }
--- a/Show More
+++ b/Show More