patch bump

better wording
fix
2025-08-25 20:21:06 +10:00 · 2025-07-12 18:46:57 +10:00 · 2025-07-12 11:19:26 +10:00 · 2025-07-07 15:43:29 +10:00 · 2025-07-04 19:37:49 +10:00 · 2025-06-30 12:42:43 +10:00
133 changed files with 2794 additions and 2785 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)

 ---

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

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

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

 ---
@@ -0,0 +1,97 @@
+# Try to find suitable version of SuperLU
+# Once done, this will define
+#
+#  SuperLU_FOUND        - system has SuperLU
+#  SuperLU_LIBRARY      - Link this to use SuperLU
+#  SuperLU_INCLUDE_DIR  - directory of SuperLU headers
+
+find_path(SuperLU_INCLUDE_DIR slu_ddefs.h
+  PATHS ${CMAKE_SYSTEM_INCLUDE_PATH} /usr/include /usr/local/include /opt/local/include
+  PATH_SUFFIXES superlu SuperLU ""
+)
+
+find_library(SuperLU_LIBRARY
+  NAMES superlu
+  PATHS ${CMAKE_SYSTEM_LIBRARY_PATH} /usr/lib64 /usr/lib /usr/local/lib64 /usr/local/lib /opt/local/lib64 /opt/local/lib
+)
+
+set(SuperLU_FOUND NO)
+
+if(NOT SuperLU_FIND_QUIETLY)
+  message(STATUS "Looking for suitable version of SuperLU")
+endif()
+
+if(SuperLU_LIBRARY AND SuperLU_INCLUDE_DIR)
+  
+  # determine if this is a suitable version of SuperLU
+  
+  if(EXISTS "${SuperLU_INCLUDE_DIR}/slu_ddefs.h")
+    
+    file(READ "${SuperLU_INCLUDE_DIR}/slu_ddefs.h" _slu_HEADER_CONTENTS)
+    
+    string(REGEX REPLACE ".*version ([0-9]).*" "\\1" SLU_VERSION_MAJOR "${_slu_HEADER_CONTENTS}")
+    
+    if("${SLU_VERSION_MAJOR}" EQUAL "7")
+      set(SuperLU_FOUND YES)
+    endif()
+    
+    if("${SLU_VERSION_MAJOR}" EQUAL "6")
+      set(SuperLU_FOUND YES)
+    endif()
+    
+    if("${SLU_VERSION_MAJOR}" EQUAL "5")
+      set(SuperLU_FOUND YES)
+    endif()
+    
+    if("${SLU_VERSION_MAJOR}" EQUAL "4")
+      
+      # NOTE: SLU_VERSION_MAJOR = 4 is not 100% reliable, as the version string in slu_ddefs.h wasn't updated in SuperLU 5.0, 5.1, 5.2
+      # NOTE: We need to explicitly check the prototype of dgssvx() to determine whether this is version 4 or 5.
+      
+      string(FIND "${_slu_HEADER_CONTENTS}" "dgssvx" _slu_dgssvx_start)
+      
+      if("${_slu_dgssvx_start}" GREATER 0)
+        
+        string(SUBSTRING "${_slu_HEADER_CONTENTS}" ${_slu_dgssvx_start} -1 _slu_dgssvx_tmp)
+        
+        string(FIND "${_slu_dgssvx_tmp}" ";" _slu_dgssvx_len)
+        
+        string(SUBSTRING "${_slu_HEADER_CONTENTS}" ${_slu_dgssvx_start} ${_slu_dgssvx_len} _slu_dgssvx_proto)
+        
+        string(FIND "${_slu_dgssvx_proto}" "GlobalLU_t" _slu_dgssvx_check)
+        
+        if("${_slu_dgssvx_check}" EQUAL -1)
+          # in version 4, dgssvx() doesn't have a GlobalLU_t argument
+          if(NOT SuperLU_FIND_QUIETLY)
+            message(STATUS "Found SuperLU 4, which is insufficient")
+          endif()
+        else()
+          set(SuperLU_FOUND YES)
+        endif()
+      
+      endif()
+    
+    endif()
+  
+  else()
+  
+    message(STATUS "Could not read ${SuperLU_INCLUDE_DIR}/slu_ddefs.h to verify suitable version")
+  
+  endif()
+  
+endif()
+
+if(SuperLU_LIBRARY AND NOT SuperLU_INCLUDE_DIR)
+  if(NOT SuperLU_FIND_QUIETLY)
+    message(STATUS "Incomplete installation of SuperLU: no SuperLU headers found")
+  endif()
+endif()
+  
+if(NOT SuperLU_FOUND) 
+  if(NOT SuperLU_FIND_QUIETLY)
+    message(STATUS "Could not find suitable version of SuperLU")
+  endif()
+else()
+  message(STATUS "Found SuperLU library: ${SuperLU_LIBRARY}")
+endif()
+
@@ -1,90 +0,0 @@
-# Try to find SuperLU 5.x
-# Once done, this will define
-#
-#  SuperLU_FOUND        - system has SuperLU
-#  SuperLU_LIBRARY      - Link this to use SuperLU
-#  SuperLU_INCLUDE_DIR  - directory of SuperLU headers
-
-find_path(SuperLU_INCLUDE_DIR slu_ddefs.h
-  PATHS ${CMAKE_SYSTEM_INCLUDE_PATH} /usr/include /usr/local/include /opt/local/include
-  PATH_SUFFIXES superlu SuperLU ""
-)
-
-find_library(SuperLU_LIBRARY
-  NAMES superlu
-  PATHS ${CMAKE_SYSTEM_LIBRARY_PATH} /usr/lib64 /usr/lib /usr/local/lib64 /usr/local/lib /opt/local/lib64 /opt/local/lib
-)
-
-set(SuperLU_FOUND NO)
-
-if(NOT SuperLU_FIND_QUIETLY)
-  message(STATUS "Looking for SuperLU version 5")
-endif()
-
-if(SuperLU_LIBRARY AND SuperLU_INCLUDE_DIR)
-  
-  # Now determine if the version is version 5.
-  # This voodoo is required as SuperLU developers like to keep things messy.
-  
-  if(EXISTS "${SuperLU_INCLUDE_DIR}/slu_ddefs.h")
-    
-    file(READ "${SuperLU_INCLUDE_DIR}/slu_ddefs.h" _slu_HEADER_CONTENTS)
-    
-    string(REGEX REPLACE ".*version ([0-9]).*" "\\1" SLU_VERSION_MAJOR "${_slu_HEADER_CONTENTS}")
-    
-    if("${SLU_VERSION_MAJOR}" EQUAL "5")
-      set(SuperLU_FOUND YES)
-    else()
-      if("${SLU_VERSION_MAJOR}" EQUAL "4")
-        
-        # NOTE: SLU_VERSION_MAJOR = 4 is not 100% reliable, as the version string in slu_ddefs.h wasn't updated in SuperLU 5.0, 5.1, 5.2
-        # NOTE: We need to explicitly check the prototype of dgssvx() to determine whether this is version 4 or 5.
-        
-        string(FIND "${_slu_HEADER_CONTENTS}" "dgssvx" _slu_dgssvx_start)
-        
-        if("${_slu_dgssvx_start}" GREATER 0)
-          
-          string(SUBSTRING "${_slu_HEADER_CONTENTS}" ${_slu_dgssvx_start} -1 _slu_dgssvx_tmp)
-          
-          string(FIND "${_slu_dgssvx_tmp}" ";" _slu_dgssvx_len)
-          
-          string(SUBSTRING "${_slu_HEADER_CONTENTS}" ${_slu_dgssvx_start} ${_slu_dgssvx_len} _slu_dgssvx_proto)
-          
-          string(FIND "${_slu_dgssvx_proto}" "GlobalLU_t" _slu_dgssvx_check)
-          
-          if("${_slu_dgssvx_check}" EQUAL -1)
-            # in version 4, dgssvx() doesn't have a GlobalLU_t argument
-            if(NOT SuperLU_FIND_QUIETLY)
-              message(STATUS "Found SuperLU, but it doesn't appear to be version 5")
-            endif()
-          else()
-            set(SuperLU_FOUND YES)
-          endif()
-        endif()
-      else()
-        if(NOT SuperLU_FIND_QUIETLY)
-          message(STATUS "Found SuperLU version ${SLU_VERSION_MAJOR}, but version 5 is required")
-        endif()
-      endif()
-    endif()
-  else()
-    if(NOT SuperLU_FIND_QUIETLY)
-      message(STATUS "Found SuperLU and headers, but could not verify version 5")
-    endif()
-  endif()
-endif()
-
-if(SuperLU_LIBRARY AND NOT SuperLU_INCLUDE_DIR)
-  if(NOT SuperLU_FIND_QUIETLY)
-    message(STATUS "Incomplete installation of SuperLU: no SuperLU headers found")
-  endif()
-endif()
-  
-if(NOT SuperLU_FOUND) 
-  if(NOT SuperLU_FIND_QUIETLY)
-    message(STATUS "Could not find SuperLU")
-  endif()
-else()
-  message(STATUS "Found SuperLU: ${SuperLU_LIBRARY}")
-endif()
-
@@ -20,7 +20,7 @@ limitations under the License.
 <head>
  <meta content="text/html;charset=ISO-8859-1" http-equiv="Content-Type">
  <title>Armadillo: API Documentation</title>
-  <link rel="icon" href="armadillo_icon.png" type="image/png">
+  <link rel="icon" type="image/png" href="data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAABAAAAAQCAYAAAAf8/9hAAAAAXNSR0IArs4c6QAAAAlwSFlzAAALEwAACxMBAJqcGAAAARpJREFUOMtjYMABzFuiGYE4BohPA/FvIP6Phv8CMTcuzVxAvAGqeTYQT0HCd6AG3MSlmRWIdwNxLhY5MyC+ADVgLS4DyqAK2NDEBaAuiobKt+Ay4CZUgQRaeKwCYj8gTofKR+My4CNUQSSSWCYQT4ay+6HyprgMOADEs4BYFso3AOITQMwB5W+HGsCPMxCR2LxAfA6INZHE7gHxGwZCAOrvFUDcCsQqUKwFxP+A+CQxBpRgSTwwvIKQZm8gfgbEVkCsgYRToQZ049NsBtVshEUuA2pAAYpfoTQzEEcB8QsgtsdheA/UgFhkwV1AfBuIPwDxOyC2xeO6jVAD5kBd6gQS/AbEX4F4JRDLEwgbUHp4CsRboLETzkApAADMF8+Pl3d18QAAAABJRU5ErkJggg==">
  <style type="text/css">
  <!--
 body
@@ -150,7 +150,7 @@ hr.greyline
 <a class="noprint" style="display:scroll; position:fixed; bottom:5px; right:5px;" href="#top"><font size=-1>[top]</font></a>

 <a name="top"></a>
-<big><b>API Documentation for Armadillo 14.2</b></big>
+<big><b>API Documentation for Armadillo 14.4</b></big>
 <br>
 <br>
 <br>
@@ -197,8 +197,8 @@ Citations are useful for the continued development and maintenance of the librar
 <br>
 <br>
 Conrad Sanderson and Ryan Curtin.
-<br><i><a href="armadillo_joss_2016.pdf">Armadillo: a template-based C++ library for linear algebra</a></i>.
-<br>Journal of Open Source Software, Vol.&thinsp;1, No.&thinsp;2, pp.&thinsp;26, 2016.
+<br><i><a href="armadillo_iccae_2025.pdf">Armadillo: An Efficient Framework for Numerical Linear Algebra</a></i>.
+<br>International Conference on Computer and Automation Engineering, 2025.
 <br>
 <br>
 <!--
@@ -434,7 +434,7 @@ Conrad Sanderson and Ryan Curtin.
 <tr style="background-color: #F5F5F5;"><td><a href="#resize">resize</a></td><td>&nbsp;</td><td>change size while keeping elements and preserving layout</td></tr>
 <tr style="background-color: #F5F5F5;"><td><a href="#reverse">reverse</a></td><td>&nbsp;</td><td>reverse order of elements</td></tr>
 <tr><td><a href="#roots">roots</a></td><td>&nbsp;</td><td>roots of polynomial</td></tr>
-<tr><td><a href="#shift">shift</a></td><td>&nbsp;</td><td>shift elements</td></tr>
+<tr><td><a href="#shift">shift</a></td><td>&nbsp;</td><td>circular shift of elements</td></tr>
 <tr><td><a href="#shuffle">shuffle</a></td><td>&nbsp;</td><td>randomly shuffle elements</td></tr>
 <tr style="background-color: #F5F5F5;"><td><a href="#size">size</a></td><td>&nbsp;</td><td>obtain dimensions of given object</td></tr>
 <tr style="background-color: #F5F5F5;"><td><a href="#sort">sort</a></td><td>&nbsp;</td><td>sort elements</td></tr>
@@ -604,7 +604,7 @@ The root matrix class is <b>Mat&lt;</b><i>type</i><b>&gt;</b>, where <i>type</i>
 </li>
 <br>
 <li>
-For convenience the following typedefs have been defined:
+For convenience, the following typedefs are defined:
 <ul>
 <table style="text-align: left;" border="0" cellpadding="2" cellspacing="2">
  <tbody>
@@ -918,7 +918,7 @@ and inherits most of the member functions
 </li>
 <br>
 <li>
-For convenience the following typedefs have been defined:
+For convenience, the following typedefs are defined:
 <ul>
 <table style="text-align: left;" border="0" cellpadding="2" cellspacing="2">
  <tbody>
@@ -1217,7 +1217,7 @@ and inherits most of the member functions
 </li>
 <br>
 <li>
-For convenience the following typedefs have been defined:
+For convenience, the following typedefs are defined:
 <ul>
 <table style="text-align: left;" border="0" cellpadding="2" cellspacing="2">
  <tbody>
@@ -1475,7 +1475,7 @@ The root cube class is <b>Cube&lt;</b><i>type</i><b>&gt;</b>, where <i>type</i>
 </li>
 <br>
 <li>
-For convenience the following typedefs have been defined:
+For convenience, the following typedefs are defined:
 <ul>
 <table style="text-align: left;" border="0" cellpadding="2" cellspacing="2">
  <tbody>
@@ -1821,7 +1821,7 @@ The root sparse matrix class is <b>SpMat&lt;</b><i>type</i><b>&gt;</b>, where <i
 </li>
 <br>
 <li>
-For convenience the following typedefs have been defined:
+For convenience, the following typedefs are defined:
 <ul>
 <table style="text-align: left;" border="0" cellpadding="2" cellspacing="2">
  <tbody>
@@ -2266,10 +2266,8 @@ umat ZZ = (AA >= BB);
 See also:
 <ul>
 <li><a href="#approx_equal">approx_equal()</a></li>
-<li><a href="#pow">pow()</a></li>
 <li><a href="#any">any()</a></li>
 <li><a href="#all">all()</a></li>
-<li><a href="#affmul">affmul()</a></li>
 <li><a href="#accu">accu()</a></li>
 <li><a href="#as_scalar">as_scalar()</a></li>
 <li><a href="#find">find()</a></li>
@@ -2423,21 +2421,6 @@ An exception is thrown if the requested element is out of bounds.
 </tr>
 <tr>
 <td style="vertical-align: top;">
-<code>.at(i)</code>&nbsp;&nbsp;or&nbsp;&nbsp;<code>[i]</code>&nbsp;
-</td>
-<td style="vertical-align: top;"><br>
-</td>
-<td style="vertical-align: top;">
-As for <code>(i)</code>, but without a bounds check; not recommended; see the caveats below
-</td>
-</tr>
-<tr>
-<td>&nbsp;</td>
-<td>&nbsp;</td>
-<td>&nbsp;</td>
-</tr>
-<tr>
-<td style="vertical-align: top;">
 <code>(r,c)</code>
 </td>
 <td style="vertical-align: top;"><br>
@@ -2454,6 +2437,37 @@ An exception is thrown if the requested element is out of bounds.
 </tr>
 <tr>
 <td style="vertical-align: top;">
+<code>(r,c,s)</code>
+</td>
+<td style="vertical-align: top;"><br>
+</td>
+<td style="vertical-align: top;">
+For <i>cube</i> and <i>3D field</i> classes, access the element/object stored at row <i>r</i>, column <i>c</i>, and slice <i>s</i>.
+An exception is thrown if the requested element is out of bounds.
+</td>
+</tr>
+<tr>
+<td>&nbsp;</td>
+<td>&nbsp;</td>
+<td>&nbsp;</td>
+</tr>
+<tr>
+<td style="vertical-align: top;">
+<code>.at(i)</code>&nbsp;&nbsp;or&nbsp;&nbsp;<code>[i]</code>&nbsp;
+</td>
+<td style="vertical-align: top;"><br>
+</td>
+<td style="vertical-align: top;">
+As for <code>(i)</code>, but without a bounds check; not recommended; see the caveats below
+</td>
+</tr>
+<tr>
+<td>&nbsp;</td>
+<td>&nbsp;</td>
+<td>&nbsp;</td>
+</tr>
+<tr>
+<td style="vertical-align: top;">
 <code>.at(r,c)</code>
 </td>
 <td style="vertical-align: top;"><br>
@@ -2471,22 +2485,6 @@ As for <code>(r,c)</code>, but without a bounds check; not recommended; see the
 </tr>
 <tr>
 <td style="vertical-align: top;">
-<code>(r,c,s)</code>
-</td>
-<td style="vertical-align: top;"><br>
-</td>
-<td style="vertical-align: top;">
-For <i>cube</i> and <i>3D field</i> classes, access the element/object stored at row <i>r</i>, column <i>c</i>, and slice <i>s</i>.
-An exception is thrown if the requested element is out of bounds.
-</td>
-</tr>
-<tr>
-<td>&nbsp;</td>
-<td>&nbsp;</td>
-<td>&nbsp;</td>
-</tr>
-<tr>
-<td style="vertical-align: top;">
 <code>.at(r,c,s)</code>
 </td>
 <td style="vertical-align: top;"><br>
@@ -2550,7 +2548,6 @@ double y = v(0);
 <li><a href="#iterators_mat">iterators (dense matrices)</a></li>
 <li><a href="#iterators_cube">iterators (cubes)</a></li>
 <li><a href="#iterators_spmat">iterators (sparse matrices)</a></li>
-<li><a href="#config_hpp">config.hpp</a></li>
 </ul>
 </li>
 <br>
@@ -3410,7 +3407,7 @@ v.set_size(100);        // or:  vec v(100, fill::none);
      <td style="vertical-align: top;"><br>
      </td>
      <td>&nbsp;</td>
-      <td style="vertical-align: top;">(member function of <i>Mat</i> and <i>SpMat</i>)
+      <td style="vertical-align: top;">(member function of <i>Mat</i>, <i>SpMat</i>, <i>field</i>)
      </td>
    </tr>
    <tr>
@@ -3418,7 +3415,7 @@ v.set_size(100);        // or:  vec v(100, fill::none);
      <td style="vertical-align: top;"><br>
      </td>
      <td>&nbsp;</td>
-      <td style="vertical-align: top;">(member function of <i>Cube</i>)
+      <td style="vertical-align: top;">(member function of <i>Cube</i>, <i>field</i>)
      </td>
    </tr>
    <tr>
@@ -3426,7 +3423,7 @@ v.set_size(100);        // or:  vec v(100, fill::none);
      <td style="vertical-align: top;"><br>
      </td>
      <td>&nbsp;</td>
-      <td style="vertical-align: top;">(member function of <i>Mat</i>, <i>Cube</i>, <i>SpMat</i>)
+      <td style="vertical-align: top;">(member function of <i>Mat</i>, <i>Cube</i>, <i>SpMat</i>, <i>field</i>)
      </td>
    </tr>
  </tbody>
@@ -3512,7 +3509,7 @@ A.reshape(5,4);
      <td style="vertical-align: top;"><br>
      </td>
      <td>&nbsp;</td>
-      <td style="vertical-align: top;">(member function of <i>Mat</i> and <i>SpMat</i>)
+      <td style="vertical-align: top;">(member function of <i>Mat</i>, <i>SpMat</i>, <i>field</i>)
      </td>
    </tr>
    <tr>
@@ -3520,7 +3517,7 @@ A.reshape(5,4);
      <td style="vertical-align: top;"><br>
      </td>
      <td>&nbsp;</td>
-      <td style="vertical-align: top;">(member function of <i>Cube</i>)
+      <td style="vertical-align: top;">(member function of <i>Cube</i>, <i>field</i>)
      </td>
    </tr>
    <tr>
@@ -3528,7 +3525,7 @@ A.reshape(5,4);
      <td style="vertical-align: top;"><br>
      </td>
      <td>&nbsp;</td>
-      <td style="vertical-align: top;">(member function of <i>Mat</i>, <i>Col</i>, <i>Row</i>, <i>Cube</i>, <i>SpMat</i>)
+      <td style="vertical-align: top;">(member function of <i>Mat</i>, <i>Col</i>, <i>Row</i>, <i>Cube</i>, <i>SpMat</i>, <i>field</i>)
      </td>
    </tr>
  </tbody>
@@ -5856,10 +5853,10 @@ If the matrix expression appears to be singular, the output matrix is reset and
 <b>Caveats:</b>
 <ul>
 <li>
-if matrix <i>A</i> is know to be symmetric positive definite, <a href="#inv_sympd">inv_sympd()</a> is faster
+if the given matrix is know to be symmetric positive definite, <a href="#inv_sympd">inv_sympd()</a> is faster
 </li>
 <li>
-to solve a system of linear equations, such as <i>Z&nbsp;=&nbsp;inv(X)*Y</i>,
+to solve a system of linear equations, such as <i>X&nbsp;=&nbsp;A.i()&thinsp;*&thinsp;B</i>,
 <a href="#solve">solve()</a> can be faster and/or more accurate
 </li>
 </ul>
@@ -8210,18 +8207,13 @@ sp_mat A = sprandu&lt;sp_mat&gt;(100, 200, 0.1);
 <a name="toeplitz"></a>
 <b>toeplitz( A )</b>
 <br><b>toeplitz( A, B )</b>
-<br><b>circ_toeplitz( A )</b>
 <ul>
 <li>
-toeplitz(): generate a Toeplitz matrix, with the first column specified by <i>A</i>, and (optionally) the first row specified by <i>B</i>
+Generate a Toeplitz matrix, with the first column specified by vector <i>A</i>, and (optionally) the first row specified by vector <i>B</i>
 </li>
 <br>
 <li>
-circ_toeplitz(): generate a circulant Toeplitz matrix
-</li>
-<br>
-<li>
-A and B must be vectors
+If argument <i>B</i> is omitted, the first row is assumed to be the same as the first column
 </li>
 <br>
 <li>
@@ -8229,8 +8221,8 @@ Examples:
 <ul>
 <pre>
 vec A(5, fill::randu);
+
 mat X = toeplitz(A);
-mat Y = circ_toeplitz(A);
 </pre>
 </ul>
 </li>
@@ -8239,7 +8231,6 @@ mat Y = circ_toeplitz(A);
 <ul>
 <li><a href="https://mathworld.wolfram.com/ToeplitzMatrix.html">Toeplitz matrix in MathWorld</a></li>
 <li><a href="https://en.wikipedia.org/wiki/Toeplitz_matrix">Toeplitz matrix in Wikipedia</a></li>
-<li><a href="https://en.wikipedia.org/wiki/Circulant_matrix">Circulant matrix in Wikipedia</a></li>
 </ul>
 </li>
 <br>
@@ -9272,7 +9263,7 @@ vec b = diff(a);
 </li>
 <br>
 <li>
-<b>Caveat:</b> <a href="#norm">norm()</a> is more robust for calculating the norm, as it handles underflows and overflows
+<b>Caveat:</b> <a href="#norm">norm()</a> is preferred for calculating the 2-norm, as it reduces the incidence of numerical underflows and overflows
 </li>
 <br>
 <li>
@@ -10606,15 +10597,15 @@ The argument <i>p</i> is optional; by default <i>p&thinsp;=&thinsp;2</i> is used
 </li>
 <br>
 <li>
-For vector norm with <i>p&thinsp;=&thinsp;2</i> and matrix norm with <i>p&thinsp;=&thinsp;"fro"</i>, a robust algorithm is used to reduce the likelihood of underflows and overflows
+For vector norm with <i>p&thinsp;=&thinsp;2</i> and matrix norm with <i>p&thinsp;=&thinsp;"fro"</i>, a robust algorithm is used to reduce the incidence of numerical underflows and overflows
 </li>
 <br>
 <li>
 <b>Caveats:</b>
 <ul>
-<li>matrix 2-norm (spectral norm) is based on SVD, which is computationally intensive; a faster alternative is <a href="#norm2est">norm2est()</a></li>
 <li>to obtain the vector norm of each row or column of a matrix, use <a href="#vecnorm">vecnorm()</a></li>
 <li>to obtain the zero/Hamming pseudo-norm (number of non-zero elements), use this expression: <code><a href="#accu">accu</a>(X&nbsp;!=&nbsp;0)</code></li>
+<li>matrix 2-norm (spectral norm) is based on SVD, which is computationally intensive; a faster alternative is <a href="#norm2est">norm2est()</a></li>
 </ul>
 </li>
 <br>
@@ -11333,7 +11324,6 @@ See also:
 <li><a href="#polyval">polyval()</a></li>
 <li><a href="#polyfit">polyfit()</a></li>
 <li><a href="#imag_real">real()</a></li>
-<li><a href="#solve">solve()</a></li>
 <li><a href="https://en.wikipedia.org/wiki/Zero_of_a_function">zero of a function in Wikipedia</a></li>
 </ul>
 </li>
@@ -12172,7 +12162,11 @@ mat Y = unique(X);
 <br><b>vecnorm( X, p, dim )</b>
 <ul>
 <li>
-Compute the <i>p</i>-norm of each column vector <i>(dim&thinsp;=&thinsp;0)</i> or row vector <i>(dim&thinsp;=&thinsp;1)</i> of matrix <i>X</i>
+If <i>X</i> is a matrix, compute the <i>p</i>-norm of each column vector <i>(dim&thinsp;=&thinsp;0)</i> or row vector <i>(dim&thinsp;=&thinsp;1)</i>
+</li>
+<br>
+<li>
+If <i>X</i> is a vector and <i>dim</i> is not specified, the result is the same as <a href="#norm">norm()</a>
 </li>
 <br>
 <li>
@@ -12611,7 +12605,7 @@ For single argument functions, <i>B = trig_fn(A)</i>, where <i>trig_fn</i> is ap
 <i>tan</i>, <i>atan</i>, <i>tanh</i>, <i>atanh</i>
 </li>
 <li>
-<i>sinc</i>, defined as sinc(x) = sin(&pi;x)&thinsp;/&thinsp;(&pi;x) for x&thinsp;&ne;&thinsp;0, and sinc(x)&thinsp;=&thinsp;1 for x&thinsp;=&thinsp;0
+<i>sinc</i>, defined as sinc(x) = sin(&pi;&middot;x)&thinsp;/&thinsp;(&pi;&middot;x) for x&thinsp;&ne;&thinsp;0, and sinc(x)&thinsp;=&thinsp;1 for x&thinsp;=&thinsp;0
 </li>
 </ul>
 </li>
@@ -13058,7 +13052,7 @@ Inverse of general square matrix <i>A</i>
 <table>
 <tbody>
 <tr><td><code>inv_opts::no_ugly</code></td><td>&nbsp;&nbsp;&nbsp;</td><td>do not provide inverses for poorly conditioned matrices (where <i>rcond &lt; A.n_rows&thinsp;&middot;&thinsp;<a href="#constants">datum::eps</a></i>)</td></tr>
-<tr><td><code>inv_opts::allow_approx</code></td><td>&nbsp;&nbsp;&nbsp;</td><td>allow approximate inverses for rank deficient or poorly conditioned matrices</td></tr>
+<tr><td><code>inv_opts::allow_approx</code></td><td>&nbsp;&nbsp;&nbsp;</td><td>provide approximate inverses (similar to <a href="#pinv">pinv</a>) for rank deficient or poorly conditioned matrices</td></tr>
 </tbody>
 </table>
 </li>
@@ -13089,7 +13083,7 @@ If <i>A</i> is not square sized, a <i>std::logic_error</i> exception is thrown
 <li>if matrix <i>A</i> is know to be symmetric positive definite, <a href="#inv_sympd">inv_sympd()</a> is faster</li>
 <li>if matrix <i>A</i> is know to be diagonal, use <i>inv(&thinsp;diagmat(A)&thinsp;)</i></li>
 <li>if matrix <i>A</i> is know to be triangular, use <i>inv(&thinsp;trimatu(A)&thinsp;)</i> or <i>inv(&thinsp;trimatl(A)&thinsp;)</i></li>
-<li>to solve a system of linear equations, such as <i>Z&thinsp;=&thinsp;inv(X)*Y</i>, <a href="#solve">solve()</a> can be faster and/or more accurate</li>
+<li>to solve a system of linear equations, such as <i>Z&nbsp;=&nbsp;inv(X)&thinsp;*&thinsp;Y</i>, <a href="#solve">solve()</a> can be faster and/or more accurate</li>
 </ul>
 </li>
 <br>
@@ -13153,7 +13147,7 @@ Inverse of symmetric/hermitian positive definite matrix <i>A</i>
 <table>
 <tbody>
 <tr><td><code>inv_opts::no_ugly</code></td><td>&nbsp;&nbsp;&nbsp;</td><td>do not provide inverses for poorly conditioned matrices (where <i>rcond &lt; A.n_rows&thinsp;&middot;&thinsp;<a href="#constants">datum::eps</a></i>)</td></tr>
-<tr><td><code>inv_opts::allow_approx</code></td><td>&nbsp;&nbsp;&nbsp;</td><td>allow approximate inverses for rank deficient or poorly conditioned symmetric matrices</td></tr>
+<tr><td><code>inv_opts::allow_approx</code></td><td>&nbsp;&nbsp;&nbsp;</td><td>provide approximate inverses (similar to <a href="#pinv">pinv</a>) for rank deficient or poorly conditioned symmetric matrices</td></tr>
 </tbody>
 </table>
 </li>
@@ -13184,7 +13178,7 @@ If <i>A</i> is not square sized, a <i>std::logic_error</i> exception is thrown
 <br>
 <li>
 <b>Caveat:</b> 
-to solve a system of linear equations, such as <i>Z&nbsp;=&nbsp;inv(X)*Y</i>,
+to solve a system of linear equations, such as <i>Z&nbsp;=&nbsp;inv_sympd(X)&thinsp;*&thinsp;Y</i>,
 using <a href="#solve">solve()</a> can be faster and/or more accurate
 </li>
 <br>
@@ -13729,7 +13723,7 @@ See also:
 <br><b>solve( X, A, B, settings )</b>
 <ul>
 <li>Solve a <b>dense</b> system of linear equations, <i>A*X = B</i>, where <i>X</i> is unknown;
-similar functionality to the \ operator in Matlab/Octave, ie. <i>X&thinsp;=&thinsp;A&thinsp;\&thinsp;B</i>
+similar functionality to the \ operator in Matlab/Octave, ie. <i>X&nbsp;=&nbsp;A&thinsp;\&thinsp;B</i>
 </li>
 <br>
 <li><i>A</i> can be square sized (critically determined system), or non-square (under/over-determined system); <i>A</i> can be rank deficient
@@ -13841,7 +13835,6 @@ See also:
 <li><a href="#inv">inv()</a></li>
 <li><a href="#pinv">pinv()</a></li>
 <li><a href="#rcond">rcond()</a></li>
-<li><a href="#roots">roots()</a></li>
 <li><a href="#syl">syl()</a></li>
 <li><a href="#spsolve">spsolve()</a> - solve sparse system of linear equations</li>
 <li><a href="https://mathworld.wolfram.com/LinearSystemofEquations.html">linear system of equations in MathWorld</a></li>
@@ -19726,6 +19719,28 @@ Use "long&nbsp;long" instead of "int" when calling BLAS and LAPACK functions; th
  </tr>
  <tr>
    <td style="vertical-align: top;">
+<code>ARMA_SUPERLU_64BIT_INT</code>
+    </td>
+    <td style="vertical-align: top;">
+      &nbsp;
+    </td>
+    <td style="vertical-align: top;">
+Use "std::int64_t" instead of "int" when calling SuperLU functions
+    </td>
+  </tr>
+  <tr>
+    <td style="vertical-align: top;">
+      &nbsp;
+    </td>
+    <td style="vertical-align: top;">
+      &nbsp;
+    </td>
+    <td style="vertical-align: top;">
+      &nbsp;
+    </td>
+  </tr>
+  <tr>
+    <td style="vertical-align: top;">
 <code>ARMA_USE_FORTRAN_HIDDEN_ARGS</code>
    </td>
    <td style="vertical-align: top;">
@@ -19851,7 +19866,6 @@ You may also need to enable one or more of the following options:
 Use 64 bit integers for matrix and vector sizes.
 Automatically enabled when using a 64-bit platform, except when using Armadillo in the R environment (via RcppArmadillo).
 Useful if matrices/vectors capable of holding more than 4 billion elements are required.
-<!-- Your machine and compiler must have support for 64 bit integers (eg. via "long" or "long&nbsp;long"). -->
 This can also be enabled by adding <code>#define&nbsp;ARMA_64BIT_WORD</code> before each instance of <code>#include&nbsp;&lt;armadillo&gt;</code>.
 <br>See also the <code>ARMA_BLAS_LONG_LONG</code> option.
    </td>
@@ -20130,6 +20144,17 @@ List of additions and changes for each version:



+<a name="version_144"></a>
+<li>Version 14.4:
+<ul>
+<li>faster handling of <a href="#pow">pow()</a> and <a href="#misc_fns">square()</a> within <a href="#accu">accu()</a> and <a href="#sum">sum()</a> expressions</li>
+<li>faster <a href="#sort">sort()</a> and <a href="#sort_index">sort_index()</a> for complex matrices</li>
+<li>expanded the <a href="#field">field class</a> with <a href="#reshape_member">.reshape()</a> and <a href="#resize_member">.resize()</a> member functions</li>
+<li>more efficient handling of compound expressions by <a href="#sum">sum()</a>, <a href="#reshape">reshape()</a>, <a href="#trans">trans()</a></li>
+<li>better detection of vector expressions by <a href="#pow">pow()</a>, <a href="#imag_real">imag()</a>, <a href="#conj">conj()</a></li>
+</ul>
+</li>
+<br>
 <a name="version_142"></a>
 <li>Version 14.2:
 <ul>
@@ -21088,7 +21113,7 @@ the <i>uword</i> type then has a minimum width of 64 bits
 <li>Version 2.2:
 <ul>
 <li>added <a href="#svd_econ">svd_econ()</a></li>
-<li>added <a href="#toeplitz">circ_toeplitz()</a></li>
+<!-- <li>added <a href="#toeplitz">circ_toeplitz()</a></li> -->
 <li>added <a href="#is_vec">.is_colvec()</a> and <a href="#is_vec">.is_rowvec()</a></li>
 </ul>
 <br>
@@ -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:
@@ -130,8 +130,8 @@ struct Base
  arma_frown("use .index_min() instead") inline elem_type min(uword& index_of_min_val) const;
  arma_frown("use .index_max() instead") inline elem_type max(uword& index_of_max_val) const;
  
-  arma_deprecated inline elem_type min(uword& row_of_min_val, uword& col_of_min_val) const;
-  arma_deprecated inline elem_type max(uword& row_of_max_val, uword& col_of_max_val) const;
+  arma_frown("use .index_min() with ind2sub() instead") inline elem_type min(uword& row_of_min_val, uword& col_of_min_val) const;
+  arma_frown("use .index_max() with ind2sub() instead") inline elem_type max(uword& row_of_max_val, uword& col_of_max_val) const;
  
  arma_warn_unused inline uword index_min() const;
  arma_warn_unused inline uword index_max() const;
@@ -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,8 @@ class CubeToMatOp : public Base< typename T1::elem_type, CubeToMatOp<T1, op_type
  inline          CubeToMatOp(const T1& in_m, const uword in_aux_uword);
  inline         ~CubeToMatOp();
  
-  arma_aligned const T1&   m;            //!< the operand; must be derived from BaseCube
-  arma_aligned       uword aux_uword;    //!< auxiliary data, uword format
+  const T1&   m;            //!< the operand; must be derived from BaseCube
+        uword aux_uword;    //!< auxiliary data, uword format
  
  template<typename eT2>
  constexpr bool is_alias(const Mat<eT2>&) const { return false; }
@@ -390,8 +390,8 @@ class Cube : public BaseCube< eT, Cube<eT> >
  arma_frown("use .index_min() instead") inline eT min(uword& index_of_min_val) const;
  arma_frown("use .index_max() instead") inline eT max(uword& index_of_max_val) const;
  
-  arma_deprecated inline eT min(uword& row_of_min_val, uword& col_of_min_val, uword& slice_of_min_val) const;
-  arma_deprecated inline eT max(uword& row_of_max_val, uword& col_of_max_val, uword& slice_of_max_val) const;
+  arma_frown("use .index_min() with ind2sub() instead") inline eT min(uword& row_of_min_val, uword& col_of_min_val, uword& slice_of_min_val) const;
+  arma_frown("use .index_max() with ind2sub() instead") inline eT max(uword& row_of_max_val, uword& col_of_max_val, uword& slice_of_max_val) const;
  
  arma_cold inline bool save(const std::string   name, const file_type type = arma_binary) const;
  arma_cold inline bool save(const hdf5_name&    spec, const file_type type = hdf5_binary) const;
@@ -485,6 +485,9 @@ class Cube<eT>::fixed : public Cube<eT>
  {
  private:
  
+  using Cube<eT>::mat_ptrs_local;
+  using Cube<eT>::mem_local;
+  
  static constexpr uword fixed_n_elem       = fixed_n_rows * fixed_n_cols * fixed_n_slices;
  static constexpr uword fixed_n_elem_slice = fixed_n_rows * fixed_n_cols;
  
@@ -1019,7 +1019,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 +2868,14 @@ Cube<eT>::Cube(const eOpCube<T1, eop_type>& X)
  
  init_cold();
  
+  if(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 +2898,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(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 +2928,14 @@ Cube<eT>::operator+=(const eOpCube<T1, eop_type>& X)
  
  if(bad_alias)  { const Cube<eT> tmp(X); return (*this).operator+=(tmp); }
  
+  if(is_same_type<eop_type, eop_pow>::value)
+    {
+    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 +2958,14 @@ Cube<eT>::operator-=(const eOpCube<T1, eop_type>& X)
  
  if(bad_alias)  { const Cube<eT> tmp(X); return (*this).operator-=(tmp); }
  
+  if(is_same_type<eop_type, eop_pow>::value)
+    {
+    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 +2988,14 @@ Cube<eT>::operator%=(const eOpCube<T1, eop_type>& X)
  
  if(bad_alias)  { const Cube<eT> tmp(X); return (*this).operator%=(tmp); }
  
+  if(is_same_type<eop_type, eop_pow>::value)
+    {
+    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 +3018,14 @@ Cube<eT>::operator/=(const eOpCube<T1, eop_type>& X)
  
  if(bad_alias)  { const Cube<eT> tmp(X); return (*this).operator/=(tmp); }
  
+  if(is_same_type<eop_type, eop_pow>::value)
+    {
+    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 +3496,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 +3509,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 +3565,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 +3585,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 +5460,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();
@@ -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();  }
@@ -970,18 +970,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 +1008,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 +1478,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 +1490,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 +1540,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();
        }
      }
    }
@@ -544,8 +544,8 @@ class Mat : public Base< eT, Mat<eT> >
  arma_frown("use .index_min() instead") inline eT min(uword& index_of_min_val) const;
  arma_frown("use .index_max() instead") inline eT max(uword& index_of_max_val) const;
  
-  arma_deprecated inline eT min(uword& row_of_min_val, uword& col_of_min_val) const;
-  arma_deprecated inline eT max(uword& row_of_max_val, uword& col_of_max_val) const;
+  arma_frown("use .index_min() with ind2sub() instead") inline eT min(uword& row_of_min_val, uword& col_of_min_val) const;
+  arma_frown("use .index_max() with ind2sub() instead") inline eT max(uword& row_of_max_val, uword& col_of_max_val) const;
  
  
  arma_cold inline bool save(const std::string   name, const file_type type = arma_binary) const;
@@ -825,6 +825,8 @@ class Mat<eT>::fixed : public Mat<eT>
  {
  private:
  
+  using Mat<eT>::mem_local;
+  
  static constexpr uword fixed_n_elem = fixed_n_rows * fixed_n_cols;
  static constexpr bool  use_extra    = (fixed_n_elem > arma_config::mat_prealloc);
  
@@ -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
@@ -2067,7 +2067,7 @@ Mat<eT>::operator-=(const subview<eT>& X)



-//! in-place matrix mutiplication (using a submatrix on the right-hand-side)
+//! in-place matrix multiplication (using a submatrix on the right-hand-side)
 template<typename eT>
 inline
 Mat<eT>&
@@ -2082,7 +2082,7 @@ Mat<eT>::operator*=(const subview<eT>& X)



-//! in-place element-wise matrix mutiplication (using a submatrix on the right-hand-side)
+//! in-place element-wise matrix multiplication (using a submatrix on the right-hand-side)
 template<typename eT>
 inline
 Mat<eT>&
@@ -2257,7 +2257,7 @@ Mat<eT>::operator-=(const subview_cube<eT>& X)



-//! in-place matrix mutiplication (using a single-slice subcube on the right-hand-side)
+//! in-place matrix multiplication (using a single-slice subcube on the right-hand-side)
 template<typename eT>
 inline
 Mat<eT>&
@@ -2274,7 +2274,7 @@ Mat<eT>::operator*=(const subview_cube<eT>& X)



-//! in-place element-wise matrix mutiplication (using a single-slice subcube on the right-hand-side)
+//! in-place element-wise matrix multiplication (using a single-slice subcube on the right-hand-side)
 template<typename eT>
 inline
 Mat<eT>&
@@ -2383,7 +2383,7 @@ Mat<eT>::operator-=(const diagview<eT>& X)



-//! in-place matrix mutiplication (using a diagview on the right-hand-side)
+//! in-place matrix multiplication (using a diagview on the right-hand-side)
 template<typename eT>
 inline
 Mat<eT>&
@@ -2398,7 +2398,7 @@ Mat<eT>::operator*=(const diagview<eT>& X)



-//! in-place element-wise matrix mutiplication (using a diagview on the right-hand-side)
+//! in-place element-wise matrix multiplication (using a diagview on the right-hand-side)
 template<typename eT>
 inline
 Mat<eT>&
@@ -2862,6 +2862,8 @@ Mat<eT>::operator=(const SpSubview<eT>& X)
  
  if(X.n_rows == X.m.n_rows)
    {
+    arma_debug_print("access via arrays");
+    
    X.m.sync();
    
    const uword sv_col_start = X.aux_col1;
@@ -2889,6 +2891,8 @@ Mat<eT>::operator=(const SpSubview<eT>& X)
    }
  else
    {
+    arma_debug_print("access via iterators");
+    
    typename SpSubview<eT>::const_iterator it     = X.begin();
    typename SpSubview<eT>::const_iterator it_end = X.end();
    
@@ -2913,6 +2917,8 @@ Mat<eT>::operator+=(const SpSubview<eT>& X)
  
  if(X.n_rows == X.m.n_rows)
    {
+    arma_debug_print("access via arrays");
+    
    X.m.sync();
    
    const uword sv_col_start = X.aux_col1;
@@ -2940,6 +2946,8 @@ Mat<eT>::operator+=(const SpSubview<eT>& X)
    }
  else
    {
+    arma_debug_print("access via iterators");
+    
    typename SpSubview<eT>::const_iterator it     = X.begin();
    typename SpSubview<eT>::const_iterator it_end = X.end();
    
@@ -2964,6 +2972,8 @@ Mat<eT>::operator-=(const SpSubview<eT>& X)
  
  if(X.n_rows == X.m.n_rows)
    {
+    arma_debug_print("access via arrays");
+    
    X.m.sync();
    
    const uword sv_col_start = X.aux_col1;
@@ -2991,6 +3001,8 @@ Mat<eT>::operator-=(const SpSubview<eT>& X)
    }
  else
    {
+    arma_debug_print("access via iterators");
+    
    typename SpSubview<eT>::const_iterator it     = X.begin();
    typename SpSubview<eT>::const_iterator it_end = X.end();
    
@@ -5185,6 +5197,14 @@ Mat<eT>::Mat(const eOp<T1, eop_type>& X)
  
  init_cold();
  
+  if(is_same_type<eop_type, eop_pow>::value)
+    {
+    constexpr bool eT_non_int = is_non_integral<eT>::value;
+    
+    if(               X.aux == eT(2)   )  { eop_square::apply(*this, reinterpret_cast< const eOp<T1, eop_square>& >(X)); return; }
+    if(eT_non_int && (X.aux == eT(0.5)))  {   eop_sqrt::apply(*this, reinterpret_cast< const eOp<T1, eop_sqrt  >& >(X)); return; }
+    }
+  
  eop_type::apply(*this, X);
  }

@@ -5207,6 +5227,14 @@ Mat<eT>::operator=(const eOp<T1, eop_type>& X)
  
  init_warm(X.get_n_rows(), X.get_n_cols());
  
+  if(is_same_type<eop_type, eop_pow>::value)
+    {
+    constexpr bool eT_non_int = is_non_integral<eT>::value;
+    
+    if(               X.aux == eT(2)   )  { eop_square::apply(*this, reinterpret_cast< const eOp<T1, eop_square>& >(X)); return *this; }
+    if(eT_non_int && (X.aux == eT(0.5)))  {   eop_sqrt::apply(*this, reinterpret_cast< const eOp<T1, eop_sqrt  >& >(X)); return *this; }
+    }
+  
  eop_type::apply(*this, X);
  
  return *this;
@@ -5228,6 +5256,14 @@ Mat<eT>::operator+=(const eOp<T1, eop_type>& X)
  
  if(bad_alias)  { const Mat<eT> tmp(X); return (*this).operator+=(tmp); }
  
+  if(is_same_type<eop_type, eop_pow>::value)
+    {
+    constexpr bool eT_non_int = is_non_integral<eT>::value;
+    
+    if(               X.aux == eT(2)   )  { eop_square::apply_inplace_plus(*this, reinterpret_cast< const eOp<T1, eop_square>& >(X)); return *this; }
+    if(eT_non_int && (X.aux == eT(0.5)))  {   eop_sqrt::apply_inplace_plus(*this, reinterpret_cast< const eOp<T1, eop_sqrt  >& >(X)); return *this; }
+    }
+  
  eop_type::apply_inplace_plus(*this, X);
  
  return *this;
@@ -5249,6 +5285,14 @@ Mat<eT>::operator-=(const eOp<T1, eop_type>& X)
  
  if(bad_alias)  { const Mat<eT> tmp(X); return (*this).operator-=(tmp); }
  
+  if(is_same_type<eop_type, eop_pow>::value)
+    {
+    constexpr bool eT_non_int = is_non_integral<eT>::value;
+    
+    if(               X.aux == eT(2)   )  { eop_square::apply_inplace_minus(*this, reinterpret_cast< const eOp<T1, eop_square>& >(X)); return *this; }
+    if(eT_non_int && (X.aux == eT(0.5)))  {   eop_sqrt::apply_inplace_minus(*this, reinterpret_cast< const eOp<T1, eop_sqrt  >& >(X)); return *this; }
+    }
+  
  eop_type::apply_inplace_minus(*this, X);
  
  return *this;
@@ -5287,6 +5331,14 @@ Mat<eT>::operator%=(const eOp<T1, eop_type>& X)
  
  if(bad_alias)  { const Mat<eT> tmp(X); return (*this).operator%=(tmp); }
  
+  if(is_same_type<eop_type, eop_pow>::value)
+    {
+    constexpr bool eT_non_int = is_non_integral<eT>::value;
+    
+    if(               X.aux == eT(2)   )  { eop_square::apply_inplace_schur(*this, reinterpret_cast< const eOp<T1, eop_square>& >(X)); return *this; }
+    if(eT_non_int && (X.aux == eT(0.5)))  {   eop_sqrt::apply_inplace_schur(*this, reinterpret_cast< const eOp<T1, eop_sqrt  >& >(X)); return *this; }
+    }
+  
  eop_type::apply_inplace_schur(*this, X);
  
  return *this;
@@ -5308,6 +5360,14 @@ Mat<eT>::operator/=(const eOp<T1, eop_type>& X)
  
  if(bad_alias)  { const Mat<eT> tmp(X); return (*this).operator/=(tmp); }
  
+  if(is_same_type<eop_type, eop_pow>::value)
+    {
+    constexpr bool eT_non_int = is_non_integral<eT>::value;
+    
+    if(               X.aux == eT(2)   )  { eop_square::apply_inplace_div(*this, reinterpret_cast< const eOp<T1, eop_square>& >(X)); return *this; }
+    if(eT_non_int && (X.aux == eT(0.5)))  {   eop_sqrt::apply_inplace_div(*this, reinterpret_cast< const eOp<T1, eop_sqrt  >& >(X)); return *this; }
+    }
+  
  eop_type::apply_inplace_div(*this, X);
  
  return *this;
@@ -6391,7 +6451,7 @@ Mat<eT>::at_alt(const uword ii) const



-//! linear element accessor (treats the matrix as a vector); bounds checking not done when ARMA_NO_DEBUG is defined
+//! linear element accessor (treats the matrix as a vector); bounds checking not done when ARMA_DONT_CHECK_CONFORMANCE is defined
 template<typename eT>
 arma_inline
 eT&
@@ -6404,7 +6464,7 @@ Mat<eT>::operator() (const uword ii)



-//! linear element accessor (treats the matrix as a vector); bounds checking not done when ARMA_NO_DEBUG is defined
+//! linear element accessor (treats the matrix as a vector); bounds checking not done when ARMA_DONT_CHECK_CONFORMANCE is defined
 template<typename eT>
 arma_inline
 const eT&
@@ -6460,7 +6520,7 @@ Mat<eT>::at(const uword ii) const



-//! element accessor; bounds checking not done when ARMA_NO_DEBUG is defined
+//! element accessor; bounds checking not done when ARMA_DONT_CHECK_CONFORMANCE is defined
 template<typename eT>
 arma_inline
 eT&
@@ -6473,7 +6533,7 @@ Mat<eT>::operator() (const uword in_row, const uword in_col)



-//! element accessor; bounds checking not done when ARMA_NO_DEBUG is defined
+//! element accessor; bounds checking not done when ARMA_DONT_CHECK_CONFORMANCE is defined
 template<typename eT>
 arma_inline
 const eT&
@@ -6626,7 +6686,7 @@ Mat<eT>::is_colvec() const



-//! returns true if the object has the same number of non-zero rows and columnns
+//! returns true if the object has the same number of non-zero rows and columns
 template<typename eT>
 arma_inline
 bool
@@ -9565,11 +9625,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
  };


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


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


@@ -79,8 +79,8 @@ struct SpBase
  arma_frown("use .index_min() instead") inline elem_type min(uword& index_of_min_val) const;
  arma_frown("use .index_max() instead") inline elem_type max(uword& index_of_max_val) const;
  
-  arma_deprecated inline elem_type min(uword& row_of_min_val, uword& col_of_min_val) const;
-  arma_deprecated inline elem_type max(uword& row_of_max_val, uword& col_of_max_val) const;
+  arma_frown("use .index_min() with ind2sub() instead") inline elem_type min(uword& row_of_min_val, uword& col_of_min_val) const;
+  arma_frown("use .index_max() with ind2sub() instead") inline elem_type max(uword& row_of_max_val, uword& col_of_max_val) const;
  
  arma_warn_unused inline uword index_min() const;
  arma_warn_unused inline uword index_max() const;
@@ -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
@@ -5202,7 +5195,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 +5474,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 +5529,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 +5911,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 +5929,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();
@@ -42,10 +42,10 @@ class SpToDOp : public Base< typename T1::elem_type, SpToDOp<T1, op_type> >
  template<typename eT2>
  constexpr bool is_alias(const Mat<eT2>&) const { return false; }
  
-  arma_aligned const T1&       m;            //!< the operand; must be derived from SpBase
-  arma_aligned       elem_type aux;          //!< auxiliary data, using the element type as used by T1
-  arma_aligned       uword     aux_uword_a;  //!< auxiliary data, uword format
-  arma_aligned       uword     aux_uword_b;  //!< auxiliary data, uword format
+  const T1&       m;            //!< the operand; must be derived from SpBase
+        elem_type aux;          //!< auxiliary data, using the element type as used by T1
+        uword     aux_uword_a;  //!< auxiliary data, uword format
+        uword     aux_uword_b;  //!< auxiliary data, uword format
  };


@@ -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;
@@ -130,6 +130,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;
@@ -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 2
-#define ARMA_VERSION_PATCH 0
-#define ARMA_VERSION_NAME  "Smooth Caffeine"
+#define ARMA_VERSION_MINOR 4
+#define ARMA_VERSION_PATCH 4
+#define ARMA_VERSION_NAME  "Filtered Espresso"



@@ -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)
@@ -7668,6 +7668,8 @@ inline
 void_ptr
 ptr_cast(blas_int (*function)(const T*, const T*, const T*))
  {
+  // TODO: investigate replacement of union-based conversion
+  
  union converter
    {
    blas_int (*fn)(const T*, const T*, const T*);
@@ -7689,6 +7691,8 @@ inline
 void_ptr
 ptr_cast(blas_int (*function)(const std::complex<T>*, const std::complex<T>*))
  {
+  // TODO: investigate replacement of union-based conversion
+  
  union converter
    {
    blas_int (*fn)(const std::complex<T>*, const std::complex<T>*);
@@ -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
  
@@ -169,7 +169,7 @@
  // https://gcc.gnu.org/bugzilla/show_bug.cgi?id=57580
  
  #if (ARMA_GCC_VERSION < 80100)
-    #pragma message("INFO: support for GCC versions older than 8.1 is deprecated"
+    #pragma message("INFO: support for GCC versions older than 8.1 is deprecated")
  #endif
  
  #define ARMA_GOOD_COMPILER
@@ -477,10 +477,6 @@
 #undef major


-// optionally allow disabling of compile-time deprecation messages (not recommended)
-// NOTE: option 'ARMA_IGNORE_DEPRECATED_MARKER' will be removed
-// NOTE: disabling deprecation messages is counter-productive
-
 #if defined(ARMA_IGNORE_DEPRECATED_MARKER)
  #undef  arma_deprecated
  #define arma_deprecated
@@ -211,7 +211,7 @@


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

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

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

+#if (defined(ARMA_BLAS_LONG_LONG) && defined(ARMA_USE_WRAPPER))
+  #pragma message ("WARNING: use of ARMA_BLAS_LONG_LONG in conjunction with ARMA_USE_WRAPPER is not supported")
+#endif
+
+#if (defined(ARMA_BLAS_64BIT_INT) && defined(ARMA_USE_WRAPPER))
+  #pragma message ("WARNING: use of ARMA_BLAS_64BIT_INT in conjunction with ARMA_USE_WRAPPER is not supported")
+#endif
+
+#if (defined(ARMA_SUPERLU_64BIT_INT) && defined(ARMA_USE_WRAPPER))
+  #pragma message ("WARNING: use of ARMA_SUPERLU_64BIT_INT in conjunction with ARMA_USE_WRAPPER is not supported")
+#endif
+
 // for compatibility with earlier versions of Armadillo
 #if defined(ARMA_BLAS_LONG) || defined(ARMA_BLAS_LONG_LONG)
  #undef  ARMA_BLAS_64BIT_INT
@@ -366,11 +378,9 @@
  #undef ARMA_CRIPPLED_LAPACK
 #endif

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


 // if Armadillo was installed on this system via CMake and ARMA_USE_WRAPPER is not defined,
@@ -211,7 +211,7 @@


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

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

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

+#if (defined(ARMA_BLAS_LONG_LONG) && defined(ARMA_USE_WRAPPER))
+  #pragma message ("WARNING: use of ARMA_BLAS_LONG_LONG in conjunction with ARMA_USE_WRAPPER is not supported")
+#endif
+
+#if (defined(ARMA_BLAS_64BIT_INT) && defined(ARMA_USE_WRAPPER))
+  #pragma message ("WARNING: use of ARMA_BLAS_64BIT_INT in conjunction with ARMA_USE_WRAPPER is not supported")
+#endif
+
+#if (defined(ARMA_SUPERLU_64BIT_INT) && defined(ARMA_USE_WRAPPER))
+  #pragma message ("WARNING: use of ARMA_SUPERLU_64BIT_INT in conjunction with ARMA_USE_WRAPPER is not supported")
+#endif
+
 // for compatibility with earlier versions of Armadillo
 #if defined(ARMA_BLAS_LONG) || defined(ARMA_BLAS_LONG_LONG)
  #undef  ARMA_BLAS_64BIT_INT
@@ -366,11 +378,9 @@
  #undef ARMA_CRIPPLED_LAPACK
 #endif

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


 // if Armadillo was installed on this system via CMake and ARMA_USE_WRAPPER is not defined,
@@ -919,12 +919,12 @@ extern "C"
  void arma_fortran(arma_zhetri)(const char* uplo, const blas_int* n, blas_cxd* a, const blas_int* lda, const blas_int* ipiv, blas_cxd* work, blas_int* info, blas_len uplo_len) ARMA_NOEXCEPT;
  
  // rcond of symmetric matrix using pre-computed factorisation (real)
-  void arma_fortran(arma_ssycon)(const char* uplo, const blas_int* n, float*    a, const blas_int* lda, const blas_int* ipiv, const float*  anorm, float*  rcond, float*    work, blas_int* iwork, blas_int* info, blas_len uplo_len) ARMA_NOEXCEPT;
-  void arma_fortran(arma_dsycon)(const char* uplo, const blas_int* n, double*   a, const blas_int* lda, const blas_int* ipiv, const double* anorm, double* rcond, double*   work, blas_int* iwork, blas_int* info, blas_len uplo_len) ARMA_NOEXCEPT;
+  void arma_fortran(arma_ssycon)(const char* uplo, const blas_int* n, const float*    a, const blas_int* lda, const blas_int* ipiv, const float*  anorm, float*  rcond, float*    work, blas_int* iwork, blas_int* info, blas_len uplo_len) ARMA_NOEXCEPT;
+  void arma_fortran(arma_dsycon)(const char* uplo, const blas_int* n, const double*   a, const blas_int* lda, const blas_int* ipiv, const double* anorm, double* rcond, double*   work, blas_int* iwork, blas_int* info, blas_len uplo_len) ARMA_NOEXCEPT;
  
  // rcond of hermitian matrix using pre-computed factorisation (complex)
-  void arma_fortran(arma_checon)(const char* uplo, const blas_int* n, blas_cxf* a, const blas_int* lda, const blas_int* ipiv, const float*  anorm, float*  rcond, blas_cxf* work, blas_int* info, blas_len uplo_len) ARMA_NOEXCEPT;
-  void arma_fortran(arma_zhecon)(const char* uplo, const blas_int* n, blas_cxd* a, const blas_int* lda, const blas_int* ipiv, const double* anorm, double* rcond, blas_cxd* work, blas_int* info, blas_len uplo_len) ARMA_NOEXCEPT;
+  void arma_fortran(arma_checon)(const char* uplo, const blas_int* n, const blas_cxf* a, const blas_int* lda, const blas_int* ipiv, const float*  anorm, float*  rcond, blas_cxf* work, blas_int* info, blas_len uplo_len) ARMA_NOEXCEPT;
+  void arma_fortran(arma_zhecon)(const char* uplo, const blas_int* n, const blas_cxd* a, const blas_int* lda, const blas_int* ipiv, const double* anorm, double* rcond, blas_cxd* work, blas_int* info, blas_len uplo_len) ARMA_NOEXCEPT;
  
 #else
  
@@ -1275,12 +1275,12 @@ extern "C"
  void arma_fortran(arma_zhetri)(const char* uplo, const blas_int* n, blas_cxd* a, const blas_int* lda, const blas_int* ipiv, blas_cxd* work, blas_int* info) ARMA_NOEXCEPT;
  
  // rcond of symmetric matrix using pre-computed factorisation (real)
-  void arma_fortran(arma_ssycon)(const char* uplo, const blas_int* n, float*    a, const blas_int* lda, const blas_int* ipiv, const float*  anorm,  float* rcond, float*    work, blas_int* iwork, blas_int* info) ARMA_NOEXCEPT;
-  void arma_fortran(arma_dsycon)(const char* uplo, const blas_int* n, double*   a, const blas_int* lda, const blas_int* ipiv, const double* anorm, double* rcond, double*   work, blas_int* iwork, blas_int* info) ARMA_NOEXCEPT;
+  void arma_fortran(arma_ssycon)(const char* uplo, const blas_int* n, const float*    a, const blas_int* lda, const blas_int* ipiv, const float*  anorm,  float* rcond, float*    work, blas_int* iwork, blas_int* info) ARMA_NOEXCEPT;
+  void arma_fortran(arma_dsycon)(const char* uplo, const blas_int* n, const double*   a, const blas_int* lda, const blas_int* ipiv, const double* anorm, double* rcond, double*   work, blas_int* iwork, blas_int* info) ARMA_NOEXCEPT;
  
  // rcond of hermitian matrix using pre-computed factorisation (complex)
-  void arma_fortran(arma_checon)(const char* uplo, const blas_int* n, blas_cxf* a, const blas_int* lda, const blas_int* ipiv, const float*  anorm, float*  rcond, blas_cxf* work, blas_int* info) ARMA_NOEXCEPT;
-  void arma_fortran(arma_zhecon)(const char* uplo, const blas_int* n, blas_cxd* a, const blas_int* lda, const blas_int* ipiv, const double* anorm, double* rcond, blas_cxd* work, blas_int* info) ARMA_NOEXCEPT;
+  void arma_fortran(arma_checon)(const char* uplo, const blas_int* n, const blas_cxf* a, const blas_int* lda, const blas_int* ipiv, const float*  anorm, float*  rcond, blas_cxf* work, blas_int* info) ARMA_NOEXCEPT;
+  void arma_fortran(arma_zhecon)(const char* uplo, const blas_int* n, const blas_cxd* a, const blas_int* lda, const blas_int* ipiv, const double* anorm, double* rcond, blas_cxd* work, blas_int* info) ARMA_NOEXCEPT;
  
 #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;
@@ -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;
  
@@ -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,8 +38,8 @@ class eGlue : public Base< typename T1::elem_type, eGlue<T1, T2, eglue_type> >
  static constexpr bool is_row  = (Proxy<T1>::is_row  || Proxy<T2>::is_row );
  static constexpr bool is_xvec = (Proxy<T1>::is_xvec || Proxy<T2>::is_xvec);
  
-  arma_aligned const Proxy<T1> P1;
-  arma_aligned const Proxy<T2> P2;
+  const Proxy<T1> P1;
+  const Proxy<T2> P2;
  
  arma_inline ~eGlue();
  arma_inline  eGlue(const T1& in_A, const T2& in_B);
@@ -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);
@@ -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,15 @@ field<oT>::init(const field<oT>& x)
  {
  arma_debug_sigprint();
  
-  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;
-    
-    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);
-        }
-      }
-    }
+  if(this == &x)  { 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 +2180,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 +2211,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,83 @@ 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>
+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(expr.aux == eT(2))
+    {
+    typedef eOp<T1,eop_square> modified_expr_type;
+    
+    return accu( reinterpret_cast< const modified_expr_type& >(expr) );
+    }
+  
+  if((expr.aux == eT(0.5)) && is_non_integral<eT>::value)
+    {
+    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 +899,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(expr.aux == eT(2))
+    {
+    typedef eOpCube<T1,eop_square> modified_expr_type;
+    
+    return accu( reinterpret_cast< const modified_expr_type& >(expr) );
+    }
+  
+  if((expr.aux == eT(0.5)) && is_non_integral<eT>::value)
+    {
+    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 +987,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)
    {
@@ -993,6 +1119,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 +1164,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;
        }
      }
@@ -1051,6 +1199,30 @@ accu(const SpOp<T1, spop_type>& expr)
  
  if(is_vectorise)  { return accu(expr.m); }
  
+  if(is_same_type<spop_type, spop_square>::yes)
+    {
+    const SpProxy<T1> P(expr.m);
+    
+    const uword N = P.get_n_nonzero();
+    
+    if(N == 0)  { return eT(0); }
+    
+    if(SpProxy<T1>::use_iterator == false)
+      {
+      return op_dot::direct_dot(N, P.get_values(), P.get_values());
+      }
+    else
+      {
+      typename SpProxy<T1>::const_iterator_type it = P.begin();
+      
+      eT val = eT(0);
+      
+      for(uword i=0; i < N; ++i)  { const eT tmp = (*it); val += (tmp*tmp); ++it; }
+      
+      return val;
+      }
+    }
+  
  const SpMat<eT> tmp = expr;
  
  return accu(tmp);
@@ -28,6 +28,9 @@ class conv_to
  {
  public:
  
+  template<typename in_eT>
+  arma_frown("use as_scalar() instead") inline static out_eT from(const in_eT& in, const typename arma_scalar_only<in_eT>::result* junk = nullptr);
+  
  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);
  
@@ -43,6 +46,26 @@ class conv_to



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


@@ -625,12 +577,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 +619,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 +661,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 +677,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 +703,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 +716,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 +742,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 +757,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 +785,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 +801,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 +844,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 +886,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 +928,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 +970,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 +1025,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);
  }


@@ -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
      }
    }
@@ -149,8 +149,8 @@ logspace(const double A, const double B, const uword N = 50u)

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

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

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

 //! kept for compatibility with old user code
 template<typename T1>
-arma_warn_unused
+arma_frown("change arma::is_finite(X) to X.is_finite()")
 inline
 bool
 is_finite(const BaseCube<typename T1::elem_type,T1>& X)
@@ -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);
+  }
+
+
+
 //! @}
@@ -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
@@ -180,6 +180,13 @@ glue_solve_gen_full::apply(Mat<eT>& actual_out, const Base<eT,T1>& A_expr, const
    const bool is_sym    = arma_config::optimise_sym && ( (refine || equilibrate || likely_sympd || force_sym || is_band || is_triu || is_tril || auxlib::crippled_lapack(A)) ? false : is_sym_expr<T1>::eval(A_expr.get_ref()) );
    const bool try_sympd = arma_config::optimise_sym && ( (          no_sympd    || is_sym       || force_sym || is_band || is_triu || is_tril || auxlib::crippled_lapack(A)) ? false : (likely_sympd ? true : sym_helper::guess_sympd(A, uword(16))) );
    
+    arma_debug_print("glue_solve_gen_full::apply(): internal flags:");
+    arma_debug_print("is_band:   ", is_band  );
+    arma_debug_print("is_triu:   ", is_triu  );
+    arma_debug_print("is_tril:   ", is_tril  );
+    arma_debug_print("is_sym:    ", is_sym   );
+    arma_debug_print("try_sympd: ", try_sympd);
+    
    if(fast)
      {
      // fast mode: solvers without refinement and without rcond estimate
@@ -31,20 +31,20 @@ glue_times_redirect2_helper<do_inv_detect>::apply(Mat<typename T1::elem_type>& o
  
  typedef typename T1::elem_type eT;
  
-  const partial_unwrap<T1> tmp1(X.A);
-  const partial_unwrap<T2> tmp2(X.B);
+  const partial_unwrap<T1> U1(X.A);
+  const partial_unwrap<T2> U2(X.B);
  
-  const typename partial_unwrap<T1>::stored_type& A = tmp1.M;
-  const typename partial_unwrap<T2>::stored_type& B = tmp2.M;
+  const typename partial_unwrap<T1>::stored_type& A = U1.M;
+  const typename partial_unwrap<T2>::stored_type& B = U2.M;
  
  constexpr bool use_alpha = partial_unwrap<T1>::do_times || partial_unwrap<T2>::do_times;
-  const     eT       alpha = use_alpha ? (tmp1.get_val() * tmp2.get_val()) : eT(0);
+  const     eT       alpha = use_alpha ? (U1.get_val() * U2.get_val()) : eT(0);
  
  if( (is_cx<eT>::no) && (resolves_to_rowvector<T1>::value && resolves_to_colvector<T2>::value) )
    {
    arma_debug_print("glue_times: dot product optimisation");
    
-    arma_conform_assert_mul_size(A, B, tmp1.do_trans, tmp2.do_trans, "matrix multiplication");
+    arma_conform_assert_mul_size(A, B, U1.do_trans, U2.do_trans, "matrix multiplication");
    
    const eT val = op_dot::direct_dot(A.n_elem, A.memptr(), B.memptr());
    
@@ -55,7 +55,7 @@ glue_times_redirect2_helper<do_inv_detect>::apply(Mat<typename T1::elem_type>& o
    return;
    }
  
-  const bool alias = tmp1.is_alias(out) || tmp2.is_alias(out);
+  const bool alias = U1.is_alias(out) || U2.is_alias(out);
  
  if(alias == false)
    {
@@ -64,7 +64,7 @@ glue_times_redirect2_helper<do_inv_detect>::apply(Mat<typename T1::elem_type>& o
      eT,
      partial_unwrap<T1>::do_trans,
      partial_unwrap<T2>::do_trans,
-      (partial_unwrap<T1>::do_times || partial_unwrap<T2>::do_times)
+      use_alpha
      >
      (out, A, B, alpha);
    }
@@ -77,7 +77,7 @@ glue_times_redirect2_helper<do_inv_detect>::apply(Mat<typename T1::elem_type>& o
      eT,
      partial_unwrap<T1>::do_trans,
      partial_unwrap<T2>::do_trans,
-      (partial_unwrap<T1>::do_times || partial_unwrap<T2>::do_times)
+      use_alpha
      >
      (tmp, A, B, alpha);
    
@@ -104,6 +104,40 @@ glue_times_redirect2_helper<true>::apply(Mat<typename T1::elem_type>& out, const
    
    const strip_inv<T1> A_strip(X.A);
    
+    typedef typename strip_inv<T1>::stored_type T1_stripped;
+    
+    if( (is_cx<eT>::no) && (strip_inv<T1>::do_inv_gen) && (is_Mat<T1_stripped>::value) && (is_Mat<T2>::value) )
+      {
+      const unwrap<T1_stripped> UA(A_strip.M);
+      const unwrap<T2         > UB(X.B);
+      
+      const typename unwrap<T1_stripped>::stored_type& A = UA.M;
+      const typename unwrap<T2         >::stored_type& B = UB.M;
+      
+      const uword N = A.n_rows;
+      
+      if( (N > 0) && (N <= uword(3)) && (N == A.n_cols) && (N == B.n_rows) && (void_ptr(&out) != void_ptr(&B)) )
+        {
+        arma_debug_print("glue_times_redirect<2>::apply(): inv tiny matrix optimisation");
+        
+        Mat<eT> AA(N, N, arma_nozeros_indicator());
+        
+        arrayops::copy(AA.memptr(), A.memptr(), AA.n_elem);
+        
+        bool inv_status = false;
+        
+        if(N == 1)  { const eT a = AA[0]; AA[0] = eT(1) / a; inv_status = (a != eT(0)); }
+        if(N == 2)  { inv_status = op_inv_gen_full::apply_tiny_2x2(AA); }
+        if(N == 3)  { inv_status = op_inv_gen_full::apply_tiny_3x3(AA); }
+        
+        if(inv_status)  { glue_times::apply<eT,false,false,false>(out, AA, B, eT(0)); return; }
+        
+        arma_debug_print("glue_times_redirect<2>::apply(): inv tiny matrix optimisation failed");
+        
+        // fallthrough if optimisation failed
+        }
+      }
+    
    Mat<eT> A = A_strip.M;
    
    arma_conform_check( (A.is_square() == false), "inv(): given matrix must be square sized" );
@@ -186,18 +220,18 @@ glue_times_redirect3_helper<do_inv_detect>::apply(Mat<typename T1::elem_type>& o
  // we have exactly 3 objects
  // hence we can safely expand X as X.A.A, X.A.B and X.B
  
-  const partial_unwrap<T1> tmp1(X.A.A);
-  const partial_unwrap<T2> tmp2(X.A.B);
-  const partial_unwrap<T3> tmp3(X.B  );
+  const partial_unwrap<T1> U1(X.A.A);
+  const partial_unwrap<T2> U2(X.A.B);
+  const partial_unwrap<T3> U3(X.B  );
  
-  const typename partial_unwrap<T1>::stored_type& A = tmp1.M;
-  const typename partial_unwrap<T2>::stored_type& B = tmp2.M;
-  const typename partial_unwrap<T3>::stored_type& C = tmp3.M;
+  const typename partial_unwrap<T1>::stored_type& A = U1.M;
+  const typename partial_unwrap<T2>::stored_type& B = U2.M;
+  const typename partial_unwrap<T3>::stored_type& C = U3.M;
  
  constexpr bool use_alpha = partial_unwrap<T1>::do_times || partial_unwrap<T2>::do_times || partial_unwrap<T3>::do_times;
-  const     eT       alpha = use_alpha ? (tmp1.get_val() * tmp2.get_val() * tmp3.get_val()) : eT(0);
+  const     eT       alpha = use_alpha ? (U1.get_val() * U2.get_val() * U3.get_val()) : eT(0);
  
-  const bool alias = tmp1.is_alias(out) || tmp2.is_alias(out) || tmp3.is_alias(out);
+  const bool alias = U1.is_alias(out) || U2.is_alias(out) || U3.is_alias(out);
  
  if(alias == false)
    {
@@ -207,7 +241,7 @@ glue_times_redirect3_helper<do_inv_detect>::apply(Mat<typename T1::elem_type>& o
      partial_unwrap<T1>::do_trans,
      partial_unwrap<T2>::do_trans,
      partial_unwrap<T3>::do_trans,
-      (partial_unwrap<T1>::do_times || partial_unwrap<T2>::do_times || partial_unwrap<T3>::do_times)
+      use_alpha
      >
      (out, A, B, C, alpha);
    }
@@ -221,7 +255,7 @@ glue_times_redirect3_helper<do_inv_detect>::apply(Mat<typename T1::elem_type>& o
      partial_unwrap<T1>::do_trans,
      partial_unwrap<T2>::do_trans,
      partial_unwrap<T3>::do_trans,
-      (partial_unwrap<T1>::do_times || partial_unwrap<T2>::do_times || partial_unwrap<T3>::do_times)
+      use_alpha
      >
      (tmp, A, B, C, alpha);
    
@@ -252,14 +286,14 @@ glue_times_redirect3_helper<true>::apply(Mat<typename T1::elem_type>& out, const
    
    arma_conform_check( (A.is_square() == false), "inv(): given matrix must be square sized" );
    
-    const partial_unwrap<T2> tmp2(X.A.B);
-    const partial_unwrap<T3> tmp3(X.B  );
+    const partial_unwrap<T2> U2(X.A.B);
+    const partial_unwrap<T3> U3(X.B  );
    
-    const typename partial_unwrap<T2>::stored_type& B = tmp2.M;
-    const typename partial_unwrap<T3>::stored_type& C = tmp3.M;
+    const typename partial_unwrap<T2>::stored_type& B = U2.M;
+    const typename partial_unwrap<T3>::stored_type& C = U3.M;
    
    constexpr bool use_alpha = partial_unwrap<T2>::do_times || partial_unwrap<T3>::do_times;
-    const     eT       alpha = use_alpha ? (tmp2.get_val() * tmp3.get_val()) : eT(0);
+    const     eT       alpha = use_alpha ? (U2.get_val() * U3.get_val()) : eT(0);
    
    Mat<eT> BC;
    
@@ -268,7 +302,7 @@ glue_times_redirect3_helper<true>::apply(Mat<typename T1::elem_type>& out, const
      eT,
      partial_unwrap<T2>::do_trans,
      partial_unwrap<T3>::do_trans,
-      (partial_unwrap<T2>::do_times || partial_unwrap<T3>::do_times)
+      use_alpha
      >
      (BC, B, C, alpha);
    
@@ -306,8 +340,8 @@ glue_times_redirect3_helper<true>::apply(Mat<typename T1::elem_type>& out, const
    
    arma_conform_check( (B.is_square() == false), "inv(): given matrix must be square sized" );
    
-    const unwrap<T3> C_tmp(X.B);
-    const Mat<eT>& C = C_tmp.M;
+    const quasi_unwrap<T3> U3(X.B);
+    const Mat<eT>& C =     U3.M;
    
    arma_conform_assert_mul_size(B, C, "matrix multiplication");
    
@@ -330,21 +364,25 @@ glue_times_redirect3_helper<true>::apply(Mat<typename T1::elem_type>& out, const
      return;
      }
    
-    const partial_unwrap_check<T1> tmp1(X.A.A, out);
+    const partial_unwrap<T1> U1(X.A.A);
    
-    const typename partial_unwrap_check<T1>::stored_type& A = tmp1.M;
+    const typename partial_unwrap<T1>::stored_type& A = U1.M;
    
-    constexpr bool use_alpha = partial_unwrap_check<T1>::do_times;
-    const     eT       alpha = use_alpha ? tmp1.get_val() : eT(0);
+    constexpr bool use_alpha = partial_unwrap<T1>::do_times;
+    const     eT       alpha = use_alpha ? U1.get_val() : eT(0);
    
-    glue_times::apply
-      <
-      eT,
-      partial_unwrap_check<T1>::do_trans,
-      false,
-      partial_unwrap_check<T1>::do_times
-      >
-      (out, A, solve_result, alpha);
+    if(U1.is_alias(out))
+      {
+      Mat<eT> tmp;
+      
+      glue_times::apply<eT, partial_unwrap<T1>::do_trans, false, use_alpha>(tmp, A, solve_result, alpha);
+      
+      out.steal_mem(tmp);
+      }
+    else
+      {
+      glue_times::apply<eT, partial_unwrap<T1>::do_trans, false, use_alpha>(out, A, solve_result, alpha);
+      }
    
    return;
    }
@@ -365,16 +403,16 @@ glue_times_redirect<N>::apply(Mat<typename T1::elem_type>& out, const Glue<T1,T2
  
  typedef typename T1::elem_type eT;
  
-  const partial_unwrap<T1> tmp1(X.A);
-  const partial_unwrap<T2> tmp2(X.B);
+  const partial_unwrap<T1> U1(X.A);
+  const partial_unwrap<T2> U2(X.B);
  
-  const typename partial_unwrap<T1>::stored_type& A = tmp1.M;
-  const typename partial_unwrap<T2>::stored_type& B = tmp2.M;
+  const typename partial_unwrap<T1>::stored_type& A = U1.M;
+  const typename partial_unwrap<T2>::stored_type& B = U2.M;
  
  constexpr bool use_alpha = partial_unwrap<T1>::do_times || partial_unwrap<T2>::do_times;
-  const     eT       alpha = use_alpha ? (tmp1.get_val() * tmp2.get_val()) : eT(0);
+  const     eT       alpha = use_alpha ? (U1.get_val() * U2.get_val()) : eT(0);
  
-  const bool alias = tmp1.is_alias(out) || tmp2.is_alias(out);
+  const bool alias = U1.is_alias(out) || U2.is_alias(out);
  
  if(alias == false)
    {
@@ -383,7 +421,7 @@ glue_times_redirect<N>::apply(Mat<typename T1::elem_type>& out, const Glue<T1,T2
      eT,
      partial_unwrap<T1>::do_trans,
      partial_unwrap<T2>::do_trans,
-      (partial_unwrap<T1>::do_times || partial_unwrap<T2>::do_times)
+      use_alpha
      >
      (out, A, B, alpha);
    }
@@ -396,7 +434,7 @@ glue_times_redirect<N>::apply(Mat<typename T1::elem_type>& out, const Glue<T1,T2
      eT,
      partial_unwrap<T1>::do_trans,
      partial_unwrap<T2>::do_trans,
-      (partial_unwrap<T1>::do_times || partial_unwrap<T2>::do_times)
+      use_alpha
      >
      (tmp, A, B, alpha);
    
@@ -446,20 +484,20 @@ glue_times_redirect<4>::apply(Mat<typename T1::elem_type>& out, const Glue< Glue
  // there is exactly 4 objects
  // hence we can safely expand X as X.A.A.A, X.A.A.B, X.A.B and X.B
  
-  const partial_unwrap<T1> tmp1(X.A.A.A);
-  const partial_unwrap<T2> tmp2(X.A.A.B);
-  const partial_unwrap<T3> tmp3(X.A.B  );
-  const partial_unwrap<T4> tmp4(X.B    );
+  const partial_unwrap<T1> U1(X.A.A.A);
+  const partial_unwrap<T2> U2(X.A.A.B);
+  const partial_unwrap<T3> U3(X.A.B  );
+  const partial_unwrap<T4> U4(X.B    );
  
-  const typename partial_unwrap<T1>::stored_type& A = tmp1.M;
-  const typename partial_unwrap<T2>::stored_type& B = tmp2.M;
-  const typename partial_unwrap<T3>::stored_type& C = tmp3.M;
-  const typename partial_unwrap<T4>::stored_type& D = tmp4.M;
+  const typename partial_unwrap<T1>::stored_type& A = U1.M;
+  const typename partial_unwrap<T2>::stored_type& B = U2.M;
+  const typename partial_unwrap<T3>::stored_type& C = U3.M;
+  const typename partial_unwrap<T4>::stored_type& D = U4.M;
  
  constexpr bool use_alpha = partial_unwrap<T1>::do_times || partial_unwrap<T2>::do_times || partial_unwrap<T3>::do_times || partial_unwrap<T4>::do_times;
-  const     eT       alpha = use_alpha ? (tmp1.get_val() * tmp2.get_val() * tmp3.get_val() * tmp4.get_val()) : eT(0);
+  const     eT       alpha = use_alpha ? (U1.get_val() * U2.get_val() * U3.get_val() * U4.get_val()) : eT(0);
  
-  const bool alias = tmp1.is_alias(out) || tmp2.is_alias(out) || tmp3.is_alias(out) || tmp4.is_alias(out);
+  const bool alias = U1.is_alias(out) || U2.is_alias(out) || U3.is_alias(out) || U4.is_alias(out);
  
  if(alias == false)
    {
@@ -470,7 +508,7 @@ glue_times_redirect<4>::apply(Mat<typename T1::elem_type>& out, const Glue< Glue
      partial_unwrap<T2>::do_trans,
      partial_unwrap<T3>::do_trans,
      partial_unwrap<T4>::do_trans,
-      (partial_unwrap<T1>::do_times || partial_unwrap<T2>::do_times || partial_unwrap<T3>::do_times || partial_unwrap<T4>::do_times)
+      use_alpha
      >
      (out, A, B, C, D, alpha);
    }
@@ -485,7 +523,7 @@ glue_times_redirect<4>::apply(Mat<typename T1::elem_type>& out, const Glue< Glue
      partial_unwrap<T2>::do_trans,
      partial_unwrap<T3>::do_trans,
      partial_unwrap<T4>::do_trans,
-      (partial_unwrap<T1>::do_times || partial_unwrap<T2>::do_times || partial_unwrap<T3>::do_times || partial_unwrap<T4>::do_times)
+      use_alpha
      >
      (tmp, A, B, C, D, alpha);
    
@@ -504,7 +542,7 @@ glue_times::apply(Mat<typename T1::elem_type>& out, const Glue<T1,T2,glue_times>
  
  constexpr uword N_mat = 1 + depth_lhs< glue_times, Glue<T1,T2,glue_times> >::num;
  
-  arma_debug_print(arma_str::format("N_mat: %u") % N_mat);
+  arma_debug_print(arma_str::format("glue_times::apply(): N_mat: %u") % N_mat);
  
  glue_times_redirect<N_mat>::apply(out, X);
  }
@@ -518,7 +556,11 @@ glue_times::apply_inplace(Mat<typename T1::elem_type>& out, const T1& X)
  {
  arma_debug_sigprint();
  
-  out = out * X;
+  typedef typename T1::elem_type eT;
+  
+  Mat<eT> tmp = out * X;
+  
+  out.steal_mem(tmp);
  }


@@ -533,9 +575,9 @@ glue_times::apply_inplace_plus(Mat<typename T1::elem_type>& out, const Glue<T1,
  typedef typename T1::elem_type            eT;
  typedef typename get_pod_type<eT>::result  T;
  
-  if( (is_outer_product<T1>::value) || (has_op_inv_any<T1>::value) || (has_op_inv_any<T2>::value) )
+  if( X.is_alias(out) || (is_outer_product<T1>::value) || (has_op_inv_any<T1>::value) || (has_op_inv_any<T2>::value) )
    {
-    // partial workaround for corner cases
+    // handle aliasing and partial workaround for corner cases
    
    const Mat<eT> tmp(X);
    
@@ -544,21 +586,21 @@ glue_times::apply_inplace_plus(Mat<typename T1::elem_type>& out, const Glue<T1,
    return;
    }
  
-  const partial_unwrap_check<T1> tmp1(X.A, out);
-  const partial_unwrap_check<T2> tmp2(X.B, out);
+  const partial_unwrap<T1> U1(X.A);
+  const partial_unwrap<T2> U2(X.B);
  
-  typedef typename partial_unwrap_check<T1>::stored_type TA;
-  typedef typename partial_unwrap_check<T2>::stored_type TB;
+  typedef typename partial_unwrap<T1>::stored_type TA;
+  typedef typename partial_unwrap<T2>::stored_type TB;
  
-  const TA& A = tmp1.M;
-  const TB& B = tmp2.M;
+  const TA& A = U1.M;
+  const TB& B = U2.M;
  
-  constexpr bool do_trans_A = partial_unwrap_check<T1>::do_trans;
-  constexpr bool do_trans_B = partial_unwrap_check<T2>::do_trans;
+  constexpr bool do_trans_A = partial_unwrap<T1>::do_trans;
+  constexpr bool do_trans_B = partial_unwrap<T2>::do_trans;
  
-  const bool use_alpha = partial_unwrap_check<T1>::do_times || partial_unwrap_check<T2>::do_times || (sign < sword(0));
+  const bool use_alpha = partial_unwrap<T1>::do_times || partial_unwrap<T2>::do_times || (sign < sword(0));
  
-  const eT       alpha = use_alpha ? ( tmp1.get_val() * tmp2.get_val() * ( (sign > sword(0)) ? eT(1) : eT(-1) ) ) : eT(0);
+  const eT       alpha = use_alpha ? ( U1.get_val() * U2.get_val() * ( (sign > sword(0)) ? eT(1) : eT(-1) ) ) : eT(0);
  
  arma_conform_assert_mul_size(A, B, do_trans_A, do_trans_B, "matrix multiplication");
  
@@ -2239,7 +2239,7 @@ gmm_diag<eT>::km_iterate(const Mat<eT>& X, const uword max_iter, const bool verb
          }
        else
          {
-          // recover by using a randomly seleced sample (last resort)
+          // recover by using a randomly selected sample (last resort)
          proposed_i = as_scalar(randi<uvec>(1, distr_param(0,X_n_cols-1)));
          }
        
@@ -937,7 +937,7 @@ gmm_full<eT>::init_constants(const bool calc_chol)
          Mat<eT>& inv_fcov =  inv_fcovs.slice(g);
    
  //const bool inv_ok = auxlib::inv(tmp_inv, fcov);
-    const bool inv_ok = auxlib::inv_sympd(tmp_inv, fcov);
+    const bool inv_ok = auxlib::inv_sympd(tmp_inv, fcov);  // using inv_sympd() instead of inv() to ensure we can do cholesky decomp via op_chol
    
    eT log_det_val  = eT(0);
    eT log_det_sign = eT(0);
@@ -2270,7 +2270,7 @@ gmm_full<eT>::km_iterate(const Mat<eT>& X, const uword max_iter, const bool verb
          }
        else
          {
-          // recover by using a randomly seleced sample (last resort)
+          // recover by using a randomly selected sample (last resort)
          proposed_i = as_scalar(randi<uvec>(1, distr_param(0,X_n_cols-1)));
          }
        
@@ -16,7 +16,7 @@
 // 
 // ------------------------------------------------------------------------
 // 
-// This file includes portions of SuperLU 5.2 software,
+// This file includes portions of SuperLU 7.0 software,
 // licensed under the following conditions.
 // 
 // Copyright (c) 2003, The Regents of the University of California, through
@@ -56,7 +56,7 @@

 #undef ARMA_SLU_HEADERS_FOUND

-// Since we need to suport float, double, cx_float and cx_double,
+// Since we need to support float, double, cx_float and cx_double,
 // as well as preserve the sanity of the user,
 // we cannot simply include all the SuperLU headers due to their messy state
 // (duplicate definitions, pollution of global namespace, bizarro defines).
@@ -64,16 +64,23 @@
 // and manually specify a few SuperLU structures and function prototypes.
 //
 // CAVEAT:
-// This code requires SuperLU version 5.2,
-// and assumes that newer 5.x versions will have no API changes.
+// This code requires SuperLU version 7.0, and assumes that newer 7.x versions have no API changes.

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

@@ -84,7 +91,7 @@ namespace arma
 namespace superlu
  {
  // Include supermatrix.h.  This gives us SuperMatrix.
-  // Put it in the slu namespace.
+  // Put it in the superlu namespace.
  // For versions of SuperLU I am familiar with, supermatrix.h does not include any other files.
  // Therefore, putting it in the superlu namespace is reasonably safe.
  // This same reasoning is true for superlu_enum_consts.h.
@@ -120,7 +127,7 @@ namespace superlu
  
  #undef ARMA_SLU_STR1
  #undef ARMA_SLU_STR2
-    
+  
  #undef ARMA_SLU_HEADER_A
  #undef ARMA_SLU_HEADER_B
  
@@ -130,7 +137,7 @@ namespace superlu
      {
      int*    panel_histo;
      double* utime;
-      float*  ops;
+      float*  ops;        // NOTE: orig definition is flops_t* ops, where flops_t = float
      int     TinyPivots;
      int     RefineSteps;
      int     expansions;
@@ -174,16 +181,16 @@ namespace superlu
    
    typedef struct e_node
      {
-      int   size;
+      int_t size;
      void* mem;
      } ExpHeader;
    
    typedef struct
      {
-      int   size;
-      int   used;
-      int   top1;
-      int   top2;
+      int_t size;
+      int_t used;
+      int_t top1;
+      int_t top2;
      void* array;
      } LU_stack_t;
    
@@ -191,16 +198,16 @@ namespace superlu
      {
      int*       xsup;
      int*       supno;   
-      int*       lsub;
-      int*       xlsub;
+      int_t*     lsub;
+      int_t*     xlsub;
      void*      lusup;
-      int*       xlusup;
+      int_t*     xlusup;
      void*      ucol;
-      int*       usub;
-      int*       xusub;
-      int        nzlmax;
-      int        nzumax;
-      int        nzlumax;
+      int_t*     usub;
+      int_t*     xusub;
+      int_t      nzlmax;
+      int_t      nzumax;
+      int_t      nzlumax;
      int        n;
      LU_space_t MemModel;
      int        num_expansions;
@@ -283,23 +290,23 @@ namespace superlu
    {
    int*    panel_histo;
    double* utime;
-    float*  ops;
+    float*  ops;        // NOTE: orig definition is flops_t* ops, where flops_t = float
    int     TinyPivots;
    int     RefineSteps;
    int     expansions;
    } SuperLUStat_t;
  
-  typedef enum {NO, YES}                                          yes_no_t;
+  typedef enum {NO, YES}                                                yes_no_t;
  typedef enum {DOFACT, SamePattern, SamePattern_SameRowPerm, FACTORED} fact_t;
-  typedef enum {NOROWPERM, LargeDiag, MY_PERMR}                   rowperm_t;
+  typedef enum {NOROWPERM, LargeDiag_MC64, LargeDiag_HWPM, MY_PERMR}    rowperm_t;
  typedef enum {NATURAL, MMD_ATA, MMD_AT_PLUS_A, COLAMD,
-                METIS_AT_PLUS_A, PARMETIS, ZOLTAN, MY_PERMC}      colperm_t;
-  typedef enum {NOTRANS, TRANS, CONJ}                             trans_t;
-  typedef enum {NOREFINE, SLU_SINGLE=1, SLU_DOUBLE, SLU_EXTRA}    IterRefine_t;
-  typedef enum {SYSTEM, USER}                                     LU_space_t;
-  typedef enum {ONE_NORM, TWO_NORM, INF_NORM}                     norm_t;
-  typedef enum {SILU, SMILU_1, SMILU_2, SMILU_3}                  milu_t;
-
+                METIS_AT_PLUS_A, PARMETIS, METIS_ATA, ZOLTAN, MY_PERMC} colperm_t;
+  typedef enum {NOTRANS, TRANS, CONJ}                                   trans_t;
+  typedef enum {NOREFINE, SLU_SINGLE=1, SLU_DOUBLE, SLU_EXTRA}          IterRefine_t;
+  typedef enum {SYSTEM, USER}                                           LU_space_t;
+  typedef enum {ONE_NORM, TWO_NORM, INF_NORM}			        norm_t;
+  typedef enum {SILU, SMILU_1, SMILU_2, SMILU_3}			milu_t;
+  
  typedef struct
    {
    fact_t        Fact;
@@ -352,16 +359,16 @@ namespace superlu
  
  typedef struct e_node
    {
-    int   size;
+    int_t size;
    void* mem;
    } ExpHeader;
  
  typedef struct
    {
-    int   size;
-    int   used;
-    int   top1;
-    int   top2;
+    int_t size;
+    int_t used;
+    int_t top1;
+    int_t top2;
    void* array;
    } LU_stack_t;
  
@@ -369,16 +376,16 @@ namespace superlu
    {
    int*       xsup;
    int*       supno;   
-    int*       lsub;
-    int*       xlsub;
+    int_t*     lsub;
+    int_t*     xlsub;
    void*      lusup;
-    int*       xlusup;
+    int_t*     xlusup;
    void*      ucol;
-    int*       usub;
-    int*       xusub;
-    int        nzlmax;
-    int        nzumax;
-    int        nzlumax;
+    int_t*     usub;
+    int_t*     xusub;
+    int_t      nzlmax;
+    int_t      nzumax;
+    int_t      nzlumax;
    int        n;
    LU_space_t MemModel;
    int        num_expansions;
@@ -90,7 +90,7 @@ memory::acquire(const uword n_elem)
  #else
    {
    //return ( new(std::nothrow) eT[n_elem] );
-    out_memptr = (eT *) malloc(sizeof(eT)*n_elem);
+    out_memptr = (eT *) std::malloc(sizeof(eT)*n_elem);
    }
  #endif
  
@@ -124,7 +124,7 @@ memory::release(eT* mem)
    }
  #elif defined(ARMA_HAVE_POSIX_MEMALIGN)
    {
-    free( (void *)(mem) );
+    std::free( (void *)(mem) );
    }
  #elif defined(_MSC_VER)
    {
@@ -134,7 +134,7 @@ memory::release(eT* mem)
  #else
    {
    //delete [] mem;
-    free( (void *)(mem) );
+    std::free( (void *)(mem) );
    }
  #endif
  
@@ -33,9 +33,9 @@ class mtGlueCube : public BaseCube< out_eT, mtGlueCube<out_eT, T1, T2, glue_type
  arma_inline  mtGlueCube(const T1& in_A, const T2& in_B, const uword in_aux_uword);
  arma_inline ~mtGlueCube();
  
-  arma_aligned const T1&   A;         //!< first operand;  must be derived from BaseCube
-  arma_aligned const T2&   B;         //!< second operand; must be derived from BaseCube
-  arma_aligned       uword aux_uword; //!< storage of auxiliary data, uword format
+  const T1&   A;         //!< first operand;  must be derived from BaseCube
+  const T2&   B;         //!< second operand; must be derived from BaseCube
+        uword aux_uword; //!< storage of auxiliary data, uword format
  };


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


@@ -45,14 +45,12 @@ class mtOpCube : public BaseCube< out_eT, mtOpCube<out_eT, T1, op_type> >
  
  inline         ~mtOpCube();
  
-  
-  arma_aligned const T1&    m;            //!< the operand; must be derived from BaseCube
-  arma_aligned       in_eT  aux;          //!< auxiliary data, using the element type as used by T1
-  arma_aligned       out_eT aux_out_eT;   //!< auxiliary data, using the element type as specified by the out_eT template parameter
-  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
+        in_eT  aux;          //!< auxiliary data, using the element type as used by T1
+        out_eT aux_out_eT;   //!< auxiliary data, using the element type as specified by the out_eT template parameter
+        uword  aux_uword_a;  //!< auxiliary data, uword format
+        uword  aux_uword_b;  //!< auxiliary data, uword format
+        uword  aux_uword_c;  //!< auxiliary data, uword format
  };


@@ -51,13 +51,11 @@ class mtOp : public Base< out_eT, mtOp<out_eT, T1, op_type> >
  template<typename eT2>
  inline bool is_alias(const Mat<eT2>& X) const;
  
-  
-  arma_aligned const T1&    m;            //!< the operand; must be derived from Base
-  arma_aligned       in_eT  aux;          //!< auxiliary data, using the element type as used by T1
-  arma_aligned       out_eT aux_out_eT;   //!< auxiliary data, using the element type as specified by the out_eT template parameter
-  arma_aligned       uword  aux_uword_a;  //!< auxiliary data, uword format
-  arma_aligned       uword  aux_uword_b;  //!< auxiliary data, uword format
-  
+  const T1&    m;            //!< the operand; must be derived from Base
+        in_eT  aux;          //!< auxiliary data, using the element type as used by T1
+        out_eT aux_out_eT;   //!< auxiliary data, using the element type as specified by the out_eT template parameter
+        uword  aux_uword_a;  //!< auxiliary data, uword format
+        uword  aux_uword_b;  //!< auxiliary data, uword format
  };


@@ -44,11 +44,11 @@ class mtSpOp : public SpBase< out_eT, mtSpOp<out_eT, T1, spop_type> >
  template<typename eT2>
  arma_inline bool is_alias(const SpMat<eT2>& X) const;
  
-  arma_aligned const T1&    m;            //!< the operand; must be derived from SpBase
-  arma_aligned       in_eT  aux;          //!< auxiliary data, using the element type as used by T1
-  arma_aligned       out_eT aux_out_eT;   //!< auxiliary data, using the element type as specified by the out_eT template parameter
-  arma_aligned       uword  aux_uword_a;
-  arma_aligned       uword  aux_uword_b;
+  const T1&    m;            //!< the operand; must be derived from SpBase
+        in_eT  aux;          //!< auxiliary data, using the element type as used by T1
+        out_eT aux_out_eT;   //!< auxiliary data, using the element type as specified by the out_eT template parameter
+        uword  aux_uword_a;
+        uword  aux_uword_b;
  };


@@ -52,9 +52,9 @@ class mtSpReduceOp : public SpBase< out_eT, mtSpReduceOp<out_eT, T1, op_type> >
  inline          mtSpReduceOp(const T1& in_m, const uword in_aux_uword_a, const uword in_aux_uword_b);
  inline         ~mtSpReduceOp();
  
-  arma_aligned const T1&   m;            //!< the operand; must be derived from SpBase
-  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
+        uword aux_uword_a;  //!< auxiliary data, uword format
+        uword aux_uword_b;  //!< auxiliary data, uword format
  };


@@ -22,7 +22,7 @@


 //! \brief
-//! Matrix multplication where the matrices have differing element types.
+//! Matrix multiplication where the matrices have differing element types.
 //! Uses caching for speedup.
 //! Matrix 'C' is assumed to have been set to the correct size (ie. taking into account transposes)

@@ -236,7 +236,7 @@ class gemm_mixed_large


 //! \brief
-//! Matrix multplication where the matrices have differing element types.
+//! Matrix multiplication where the matrices have differing element types.

 template<const bool do_trans_A=false, const bool do_trans_B=false, const bool use_alpha=false, const bool use_beta=false>
 class gemm_mixed
@@ -68,9 +68,11 @@ class op_all
  static inline bool all_vec(T1& X);
  
  
-  template<typename T1>
-  static inline void apply_helper(Mat<uword>& out, const Proxy<T1>& P, const uword dim);
+  template<typename eT>
+  static inline void apply_mat_noalias(Mat<uword>& out, const Mat<eT>& X, const uword dim);
  
+  template<typename T1>
+  static inline void apply_proxy_noalias(Mat<uword>& out, const Proxy<T1>& P, const uword dim);
  
  template<typename T1>
  static inline void apply(Mat<uword>& out, const mtOp<uword, T1, op_all>& X);
@@ -277,17 +277,15 @@ op_all::all_vec(T1& X)



-template<typename T1>
+template<typename eT>
 inline
 void
-op_all::apply_helper(Mat<uword>& out, const Proxy<T1>& P, const uword dim)
+op_all::apply_mat_noalias(Mat<uword>& out, const Mat<eT>& X, const uword dim)
  {
  arma_debug_sigprint();
  
-  const uword n_rows = P.get_n_rows();
-  const uword n_cols = P.get_n_cols();
-  
-  typedef typename Proxy<T1>::elem_type eT;
+  const uword n_rows = X.n_rows;
+  const uword n_cols = X.n_cols;
  
  if(dim == 0)  // traverse rows (ie. process each column)
    {
@@ -297,37 +295,18 @@ op_all::apply_helper(Mat<uword>& out, const Proxy<T1>& P, const uword dim)
    
    uword* out_mem = out.memptr();
    
-    if(is_Mat<typename Proxy<T1>::stored_type>::value)
+    for(uword col=0; col < n_cols; ++col)
      {
-      const unwrap<typename Proxy<T1>::stored_type> U(P.Q);
+      const eT* colmem = X.colptr(col);
      
-      for(uword col=0; col < n_cols; ++col)
+      uword count = 0;
+      
+      for(uword row=0; row < n_rows; ++row)
        {
-        const eT* colmem = U.M.colptr(col);
-        
-        uword count = 0;
-        
-        for(uword row=0; row < n_rows; ++row)
-          {
-          count += (colmem[row] != eT(0)) ? uword(1) : uword(0);
-          }
-        
-        out_mem[col] = (n_rows == count) ? uword(1) : uword(0); 
-        }
-      }
-    else
-      {
-      for(uword col=0; col < n_cols; ++col)
-        {
-        uword count = 0;
-        
-        for(uword row=0; row < n_rows; ++row)
-          {
-          if(P.at(row,col) != eT(0))  { ++count; }
-          }
-        
-        out_mem[col] = (n_rows == count) ? uword(1) : uword(0); 
+        count += (colmem[row] != eT(0)) ? uword(1) : uword(0);
        }
+      
+      out_mem[col] = (n_rows == count) ? uword(1) : uword(0); 
      }
    }
  else
@@ -338,31 +317,15 @@ op_all::apply_helper(Mat<uword>& out, const Proxy<T1>& P, const uword dim)
    
    // internal dual use of 'out': keep the counts for each row
    
-    if(is_Mat<typename Proxy<T1>::stored_type>::value)
+    for(uword col=0; col < n_cols; ++col)
      {
-      const unwrap<typename Proxy<T1>::stored_type> U(P.Q);
+      const eT* colmem = X.colptr(col);
      
-      for(uword col=0; col < n_cols; ++col)
+      for(uword row=0; row < n_rows; ++row)
        {
-        const eT* colmem = U.M.colptr(col);
-        
-        for(uword row=0; row < n_rows; ++row)
-          {
-          out_mem[row] += (colmem[row] != eT(0)) ? uword(1) : uword(0);
-          }
+        out_mem[row] += (colmem[row] != eT(0)) ? uword(1) : uword(0);
        }
      }
-    else
-      {
-      for(uword col=0; col < n_cols; ++col)
-        {
-        for(uword row=0; row < n_rows; ++row)
-          {
-          if(P.at(row,col) != eT(0))  { ++out_mem[row]; }
-          }
-        }
-      }
-    
    
    // see what the counts tell us
    
@@ -370,7 +333,63 @@ op_all::apply_helper(Mat<uword>& out, const Proxy<T1>& P, const uword dim)
      {
      out_mem[row] = (n_cols == out_mem[row]) ? uword(1) : uword(0);
      }
+    }
+  }
+
+
+
+template<typename T1>
+inline
+void
+op_all::apply_proxy_noalias(Mat<uword>& out, const Proxy<T1>& P, const uword dim)
+  {
+  arma_debug_sigprint();
+  
+  typedef typename Proxy<T1>::elem_type eT;
+  
+  const uword n_rows = P.get_n_rows();
+  const uword n_cols = P.get_n_cols();
+  
+  if(dim == 0)  // traverse rows (ie. process each column)
+    {
+    out.zeros(1, n_cols);
    
+    if(out.n_elem == 0)  { return; }
+    
+    uword* out_mem = out.memptr();
+    
+    for(uword col=0; col < n_cols; ++col)
+      {
+      uword count = 0;
+      
+      for(uword row=0; row < n_rows; ++row)
+        {
+        if(P.at(row,col) != eT(0))  { ++count; }
+        }
+      
+      out_mem[col] = (n_rows == count) ? uword(1) : uword(0); 
+      }
+    }
+  else
+    {
+    out.zeros(n_rows, 1);
+    
+    uword* out_mem = out.memptr();
+    
+    // internal dual use of 'out': keep the counts for each row
+    
+    for(uword col=0; col < n_cols; ++col)
+    for(uword row=0; row < n_rows; ++row)
+      {
+      if(P.at(row,col) != eT(0))  { ++out_mem[row]; }
+      }
+    
+    // see what the counts tell us
+    
+    for(uword row=0; row < n_rows; ++row)
+      {
+      out_mem[row] = (n_cols == out_mem[row]) ? uword(1) : uword(0);
+      }
    }
  }

@@ -385,19 +404,39 @@ op_all::apply(Mat<uword>& out, const mtOp<uword, T1, op_all>& X)
  
  const uword dim = X.aux_uword_a;
  
-  const Proxy<T1> P(X.m);
-  
-  if(P.is_alias(out) == false)
+  if( (is_Mat<T1>::value) || (is_Mat<typename Proxy<T1>::stored_type>::value) || (arma_config::openmp && Proxy<T1>::use_mp) )
    {
-    op_all::apply_helper(out, P, dim);
+    const quasi_unwrap<T1> U(X.m);
+    
+    if(U.is_alias(out) == false)
+      {
+      op_all::apply_mat_noalias(out, U.M, dim);
+      }
+    else
+      {
+      Mat<uword> tmp;
+      
+      op_all::apply_mat_noalias(tmp, U.M, dim);
+      
+      out.steal_mem(tmp);
+      }
    }
  else
    {
-    Mat<uword> out2;
+    const Proxy<T1> P(X.m);
    
-    op_all::apply_helper(out2, P, dim);
-    
-    out.steal_mem(out2);
+    if(P.is_alias(out) == false)
+      {
+      op_all::apply_proxy_noalias(out, P, dim);
+      }
+    else
+      {
+      Mat<uword> tmp;
+      
+      op_all::apply_proxy_noalias(tmp, P, dim);
+      
+      out.steal_mem(tmp);
+      }
    }
  }

@@ -68,9 +68,11 @@ class op_any
  static inline bool any_vec(T1& X);
  
  
-  template<typename T1>
-  static inline void apply_helper(Mat<uword>& out, const Proxy<T1>& P, const uword dim);
+  template<typename eT>
+  static inline void apply_mat_noalias(Mat<uword>& out, const Mat<eT>& X, const uword dim);
  
+  template<typename T1>
+  static inline void apply_proxy_noalias(Mat<uword>& out, const Proxy<T1>& P, const uword dim);
  
  template<typename T1>
  static inline void apply(Mat<uword>& out, const mtOp<uword, T1, op_any>& X);
@@ -269,17 +269,15 @@ op_any::any_vec(T1& X)



-template<typename T1>
+template<typename eT>
 inline
 void
-op_any::apply_helper(Mat<uword>& out, const Proxy<T1>& P, const uword dim)
+op_any::apply_mat_noalias(Mat<uword>& out, const Mat<eT>& X, const uword dim)
  {
  arma_debug_sigprint();
  
-  const uword n_rows = P.get_n_rows();
-  const uword n_cols = P.get_n_cols();
-  
-  typedef typename Proxy<T1>::elem_type eT;
+  const uword n_rows = X.n_rows;
+  const uword n_cols = X.n_cols;
  
  if(dim == 0)  // traverse rows (ie. process each column)
    {
@@ -287,28 +285,13 @@ op_any::apply_helper(Mat<uword>& out, const Proxy<T1>& P, const uword dim)
    
    uword* out_mem = out.memptr();
    
-    if(is_Mat<typename Proxy<T1>::stored_type>::value)
+    for(uword col=0; col < n_cols; ++col)
      {
-      const unwrap<typename Proxy<T1>::stored_type> U(P.Q);
+      const eT* colmem = X.colptr(col);
      
-      for(uword col=0; col < n_cols; ++col)
+      for(uword row=0; row < n_rows; ++row)
        {
-        const eT* colmem = U.M.colptr(col);
-        
-        for(uword row=0; row < n_rows; ++row)
-          {
-          if(colmem[row] != eT(0))  { out_mem[col] = uword(1); break; }
-          }
-        }
-      }
-    else
-      {
-      for(uword col=0; col < n_cols; ++col)
-        {
-        for(uword row=0; row < n_rows; ++row)
-          {
-          if(P.at(row,col) != eT(0))  { out_mem[col] = uword(1); break; }
-          }
+        if(colmem[row] != eT(0))  { out_mem[col] = uword(1); break; }
        }
      }
    }
@@ -318,28 +301,57 @@ op_any::apply_helper(Mat<uword>& out, const Proxy<T1>& P, const uword dim)
    
    uword* out_mem = out.memptr();
    
-    if(is_Mat<typename Proxy<T1>::stored_type>::value)
+    for(uword col=0; col < n_cols; ++col)
      {
-      const unwrap<typename Proxy<T1>::stored_type> U(P.Q);
+      const eT* colmem = X.colptr(col);
      
-      for(uword col=0; col < n_cols; ++col)
+      for(uword row=0; row < n_rows; ++row)
        {
-        const eT* colmem = U.M.colptr(col);
-        
-        for(uword row=0; row < n_rows; ++row)
-          {
-          if(colmem[row] != eT(0))  { out_mem[row] = uword(1); }
-          }
+        if(colmem[row] != eT(0))  { out_mem[row] = uword(1); }
        }
      }
-    else
+    }
+  }
+
+
+
+template<typename T1>
+inline
+void
+op_any::apply_proxy_noalias(Mat<uword>& out, const Proxy<T1>& P, const uword dim)
+  {
+  arma_debug_sigprint();
+  
+  typedef typename Proxy<T1>::elem_type eT;
+  
+  const uword n_rows = P.get_n_rows();
+  const uword n_cols = P.get_n_cols();
+  
+  if(dim == 0)  // traverse rows (ie. process each column)
+    {
+    out.zeros(1, n_cols);
+    
+    uword* out_mem = out.memptr();
+    
+    for(uword col=0; col < n_cols; ++col)
      {
-      for(uword col=0; col < n_cols; ++col)
+      for(uword row=0; row < n_rows; ++row)
        {
-        for(uword row=0; row < n_rows; ++row)
-          {
-          if(P.at(row,col) != eT(0))  { out_mem[row] = uword(1); }
-          }
+        if(P.at(row,col) != eT(0))  { out_mem[col] = uword(1); break; }
+        }
+      }
+    }
+  else
+    {
+    out.zeros(n_rows, 1);
+    
+    uword* out_mem = out.memptr();
+    
+    for(uword col=0; col < n_cols; ++col)
+      {
+      for(uword row=0; row < n_rows; ++row)
+        {
+        if(P.at(row,col) != eT(0))  { out_mem[row] = uword(1); }
        }
      }
    }
@@ -356,19 +368,39 @@ op_any::apply(Mat<uword>& out, const mtOp<uword, T1, op_any>& X)
  
  const uword dim = X.aux_uword_a;
  
-  const Proxy<T1> P(X.m);
-  
-  if(P.is_alias(out) == false)
+  if( (is_Mat<T1>::value) || (is_Mat<typename Proxy<T1>::stored_type>::value) || (arma_config::openmp && Proxy<T1>::use_mp) )
    {
-    op_any::apply_helper(out, P, dim);
+    const quasi_unwrap<T1> U(X.m);
+    
+    if(U.is_alias(out) == false)
+      {
+      op_any::apply_mat_noalias(out, U.M, dim);
+      }
+    else
+      {
+      Mat<uword> tmp;
+      
+      op_any::apply_mat_noalias(tmp, U.M, dim);
+      
+      out.steal_mem(tmp);
+      }
    }
  else
    {
-    Mat<uword> out2;
+    const Proxy<T1> P(X.m);
    
-    op_any::apply_helper(out2, P, dim);
-    
-    out.steal_mem(out2);
+    if(P.is_alias(out) == false)
+      {
+      op_any::apply_proxy_noalias(out, P, dim);
+      }
+    else
+      {
+      Mat<uword> tmp;
+      
+      op_any::apply_proxy_noalias(tmp, P, dim);
+      
+      out.steal_mem(tmp);
+      }
    }
  }

@@ -106,30 +106,25 @@ op_dot::direct_dot(const uword n_elem, const eT* const A, const eT* const B)
  {
  arma_debug_sigprint();
  
-  if( n_elem <= 32u )
+  if(n_elem <= 32u)  { return op_dot::direct_dot_arma(n_elem, A, B); }
+  
+  #if defined(ARMA_USE_ATLAS)
+    {
+    arma_debug_print("atlas::cblas_dot()");
+    
+    return atlas::cblas_dot(n_elem, A, B);
+    }
+  #elif defined(ARMA_USE_BLAS)
+    {
+    arma_debug_print("blas::dot()");
+    
+    return blas::dot(n_elem, A, B);
+    }
+  #else
    {
    return op_dot::direct_dot_arma(n_elem, A, B);
    }
-  else
-    {
-    #if defined(ARMA_USE_ATLAS)
-      {
-      arma_debug_print("atlas::cblas_dot()");
-      
-      return atlas::cblas_dot(n_elem, A, B);
-      }
-    #elif defined(ARMA_USE_BLAS)
-      {
-      arma_debug_print("blas::dot()");
-      
-      return blas::dot(n_elem, A, B);
-      }
-    #else
-      {
-      return op_dot::direct_dot_arma(n_elem, A, B);
-      }
-    #endif
-    }
+  #endif
  }


@@ -140,30 +135,25 @@ inline
 typename arma_cx_only<eT>::result
 op_dot::direct_dot(const uword n_elem, const eT* const A, const eT* const B)
  {
-  if( n_elem <= 16u )
+  if(n_elem <= 16u)  { return op_dot::direct_dot_arma(n_elem, A, B); }
+  
+  #if defined(ARMA_USE_ATLAS)
+    {
+    arma_debug_print("atlas::cblas_cx_dot()");
+    
+    return atlas::cblas_cx_dot(n_elem, A, B);
+    }
+  #elif defined(ARMA_USE_BLAS)
+    {
+    arma_debug_print("blas::dot()");
+    
+    return blas::dot(n_elem, A, B);
+    }
+  #else
    {
    return op_dot::direct_dot_arma(n_elem, A, B);
    }
-  else
-    {
-    #if defined(ARMA_USE_ATLAS)
-      {
-      arma_debug_print("atlas::cblas_cx_dot()");
-      
-      return atlas::cblas_cx_dot(n_elem, A, B);
-      }
-    #elif defined(ARMA_USE_BLAS)
-      {
-      arma_debug_print("blas::dot()");
-      
-      return blas::dot(n_elem, A, B);
-      }
-    #else
-      {
-      return op_dot::direct_dot_arma(n_elem, A, B);
-      }
-    #endif
-    }
+  #endif
  }


@@ -461,41 +451,36 @@ op_cdot::direct_cdot(const uword n_elem, const eT* const A, const eT* const B)
  {
  arma_debug_sigprint();
  
-  if( n_elem <= 32u )
+  if(n_elem <= 32u)  { return op_cdot::direct_cdot_arma(n_elem, A, B); }
+  
+  #if defined(ARMA_USE_BLAS)
+    {
+    arma_debug_print("blas::gemv()");
+    
+    // using gemv() workaround due to compatibility issues with cdotc() and zdotc()
+    
+    const char trans   = 'C';
+    
+    const blas_int m   = blas_int(n_elem);
+    const blas_int n   = 1;
+    //const blas_int lda = (n_elem > 0) ? blas_int(n_elem) : blas_int(1);
+    const blas_int inc = 1;
+    
+    const eT alpha     = eT(1);
+    const eT beta      = eT(0);
+    
+    eT result[2];  // paranoia: using two elements instead of one
+    
+    //blas::gemv(&trans, &m, &n, &alpha, A, &lda, B, &inc, &beta, &result[0], &inc);
+    blas::gemv(&trans, &m, &n, &alpha, A, &m, B, &inc, &beta, &result[0], &inc);
+    
+    return result[0];
+    }
+  #else
    {
    return op_cdot::direct_cdot_arma(n_elem, A, B);
    }
-  else
-    {
-    #if defined(ARMA_USE_BLAS)
-      {
-      arma_debug_print("blas::gemv()");
-      
-      // using gemv() workaround due to compatibility issues with cdotc() and zdotc()
-      
-      const char trans   = 'C';
-      
-      const blas_int m   = blas_int(n_elem);
-      const blas_int n   = 1;
-      //const blas_int lda = (n_elem > 0) ? blas_int(n_elem) : blas_int(1);
-      const blas_int inc = 1;
-      
-      const eT alpha     = eT(1);
-      const eT beta      = eT(0);
-      
-      eT result[2];  // paranoia: using two elements instead of one
-      
-      //blas::gemv(&trans, &m, &n, &alpha, A, &lda, B, &inc, &beta, &result[0], &inc);
-      blas::gemv(&trans, &m, &n, &alpha, A, &m, B, &inc, &beta, &result[0], &inc);
-      
-      return result[0];
-      }
-    #else
-      {
-      return op_cdot::direct_cdot_arma(n_elem, A, B);
-      }
-    #endif
-    }
+  #endif
  }


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

@@ -123,7 +123,7 @@ op_fft_real::apply( Mat< std::complex<typename T1::pod_type> >& out, const mtOp<
    }
  else
    {
-    // process each column seperately
+    // process each column separately
    
    out.set_size(N_user, n_cols);
    
@@ -234,7 +234,7 @@ op_fft_cx::apply_noalias(Mat<eT>& out, const Mat<eT>& X, const uword a, const uw
    }
  else
    {
-    // process each column seperately
+    // process each column separately
    
    out.set_size(N_user, n_cols);
    
@@ -306,49 +306,45 @@ op_htrans::apply_direct(Mat<typename T1::elem_type>& out, const T1& X)
  typedef typename T1::elem_type eT;
  
  // allow detection of in-place transpose
-  if(is_Mat<T1>::value || (arma_config::openmp && Proxy<T1>::use_mp))
+  if(is_Mat<T1>::value)
    {
    const unwrap<T1> U(X);
    
    op_htrans::apply_mat(out, U.M);
    }
  else
+  if((is_Mat<typename Proxy<T1>::stored_type>::value) || (arma_config::openmp && Proxy<T1>::use_mp))
    {
-    const Proxy<T1> P(X);
+    const quasi_unwrap<T1> U(X);
    
-    const bool is_alias = P.is_alias(out);
-    
-    if(is_Mat<typename Proxy<T1>::stored_type>::value)
+    if(U.is_alias(out))
      {
-      const quasi_unwrap<typename Proxy<T1>::stored_type> U(P.Q);
+      Mat<eT> tmp;
      
-      if(is_alias)
-        {
-        Mat<eT> tmp;
-        
-        op_htrans::apply_mat_noalias(tmp, U.M);
-        
-        out.steal_mem(tmp);
-        }
-      else
-        {
-        op_htrans::apply_mat_noalias(out, U.M);
-        }
+      op_htrans::apply_mat_noalias(tmp, U.M);
+      
+      out.steal_mem(tmp);
      }
    else
      {
-      if(is_alias)
-        {
-        Mat<eT> tmp;
-        
-        op_htrans::apply_proxy(tmp, P);
-        
-        out.steal_mem(tmp);
-        }
-      else
-        {
-        op_htrans::apply_proxy(out, P);
-        }
+      op_htrans::apply_mat_noalias(out, U.M);
+      }
+    }
+  else
+    {
+    const Proxy<T1> P(X);
+    
+    if(P.is_alias(out))
+      {
+      Mat<eT> tmp;
+      
+      op_htrans::apply_proxy(tmp, P);
+      
+      out.steal_mem(tmp);
+      }
+    else
+      {
+      op_htrans::apply_proxy(out, P);
      }
    }
  }
@@ -393,32 +393,19 @@ op_max::direct_max(const eT* const X, const uword n_elem, uword& index_of_max_va
  {
  arma_debug_sigprint();
  
-  eT max_val_i = priv::most_neg<eT>();
-  eT max_val_j = priv::most_neg<eT>();
+  eT    best_val   = priv::most_neg<eT>();
+  uword best_index = 0;
  
-  uword best_index_i = 0;
-  uword best_index_j = 0;
-  
-  uword i,j;
-  for(i=0, j=1; j<n_elem; i+=2, j+=2)
+  for(uword i=0; i < n_elem; ++i)
    {
-    const eT X_i = X[i];
-    const eT X_j = X[j];
+    const eT val = X[i];
    
-    if(X_i > max_val_i)  { max_val_i = X_i; best_index_i = i; }
-    if(X_j > max_val_j)  { max_val_j = X_j; best_index_j = j; }
+    if(val > best_val)  { best_val = val; best_index = i; }
    }
  
-  if(i < n_elem)
-    {
-    const eT X_i = X[i];
-    
-    if(X_i > max_val_i)  { max_val_i = X_i; best_index_i = i; }
-    }
+  index_of_max_val = best_index;
  
-  index_of_max_val = (max_val_i > max_val_j) ? best_index_i : best_index_j;
-  
-  return (max_val_i > max_val_j) ? max_val_i : max_val_j;
+  return best_val;
  }


@@ -705,6 +692,8 @@ op_max::max_with_index(const Proxy<T1>& P, uword& index_of_max_val)
    {
    arma_conform_check(true, "max(): object has no elements");
    
+    index_of_max_val = 0;
+    
    return Datum<eT>::nan;
    }
  
@@ -786,6 +775,8 @@ op_max::max_with_index(const ProxyCube<T1>& P, uword& index_of_max_val)
    {
    arma_conform_check(true, "max(): object has no elements");
    
+    index_of_max_val = 0;
+    
    return Datum<eT>::nan;
    }
  
@@ -1160,6 +1151,8 @@ op_max::max_with_index(const Proxy<T1>& P, uword& index_of_max_val)
    {
    arma_conform_check(true, "max(): object has no elements");
    
+    index_of_max_val = 0;
+    
    return Datum<eT>::nan;
    }
  
@@ -1263,6 +1256,8 @@ op_max::max_with_index(const ProxyCube<T1>& P, uword& index_of_max_val)
    {
    arma_conform_check(true, "max(): object has no elements");
    
+    index_of_max_val = 0;
+    
    return Datum<eT>::nan;
    }
  
@@ -58,7 +58,7 @@ op_mean::apply_noalias(Mat<typename T1::elem_type>& out, const Proxy<T1>& P, con
  {
  arma_debug_sigprint();
  
-  if(is_Mat<typename Proxy<T1>::stored_type>::value)
+  if((is_Mat<typename Proxy<T1>::stored_type>::value) || (arma_config::openmp && Proxy<T1>::use_mp))
    {
    op_mean::apply_noalias_unwrap(out, P, dim);
    }
@@ -245,7 +245,7 @@ op_mean::apply_noalias(Cube<typename T1::elem_type>& out, const ProxyCube<T1>& P
  {
  arma_debug_sigprint();
  
-  if(is_Cube<typename ProxyCube<T1>::stored_type>::value)
+  if((is_Cube<typename ProxyCube<T1>::stored_type>::value) || (arma_config::openmp && ProxyCube<T1>::use_mp))
    {
    op_mean::apply_noalias_unwrap(out, P, dim);
    }
@@ -41,7 +41,7 @@ op_median::apply(Mat<typename T1::elem_type>& out, const Op<T1,op_median>& expr)
    {
    Mat<eT> tmp;
    
-    op_median::apply_noalias(out, U.M, dim);
+    op_median::apply_noalias(tmp, U.M, dim);
    
    out.steal_mem(tmp);
    }
@@ -393,32 +393,19 @@ op_min::direct_min(const eT* const X, const uword n_elem, uword& index_of_min_va
  {
  arma_debug_sigprint();
  
-  eT min_val_i = priv::most_pos<eT>();
-  eT min_val_j = priv::most_pos<eT>();
+  eT    best_val   = priv::most_pos<eT>();
+  uword best_index = 0;
  
-  uword best_index_i = 0;
-  uword best_index_j = 0;
-  
-  uword i,j;
-  for(i=0, j=1; j<n_elem; i+=2, j+=2)
+  for(uword i=0; i < n_elem; ++i)
    {
-    const eT X_i = X[i];
-    const eT X_j = X[j];
+    const eT val = X[i];
    
-    if(X_i < min_val_i)  { min_val_i = X_i; best_index_i = i; }
-    if(X_j < min_val_j)  { min_val_j = X_j; best_index_j = j; }
+    if(val < best_val)  { best_val = val; best_index = i; }
    }
  
-  if(i < n_elem)
-    {
-    const eT X_i = X[i];
-    
-    if(X_i < min_val_i)  { min_val_i = X_i; best_index_i = i; }
-    }
+  index_of_min_val = best_index;
  
-  index_of_min_val = (min_val_i < min_val_j) ? best_index_i : best_index_j;
-  
-  return (min_val_i < min_val_j) ? min_val_i : min_val_j;
+  return best_val;
  }


@@ -705,6 +692,8 @@ op_min::min_with_index(const Proxy<T1>& P, uword& index_of_min_val)
    {
    arma_conform_check(true, "min(): object has no elements");
    
+    index_of_min_val = 0;
+    
    return Datum<eT>::nan;
    }
  
@@ -786,6 +775,8 @@ op_min::min_with_index(const ProxyCube<T1>& P, uword& index_of_min_val)
    {
    arma_conform_check(true, "min(): object has no elements");
    
+    index_of_min_val = 0;
+    
    return Datum<eT>::nan;
    }
  
@@ -1160,6 +1151,8 @@ op_min::min_with_index(const Proxy<T1>& P, uword& index_of_min_val)
    {
    arma_conform_check(true, "min(): object has no elements");
    
+    index_of_min_val = 0;
+    
    return Datum<eT>::nan;
    }
  
@@ -1263,6 +1256,8 @@ op_min::min_with_index(const ProxyCube<T1>& P, uword& index_of_min_val)
    {
    arma_conform_check(true, "min(): object has no elements");
    
+    index_of_min_val = 0;
+    
    return Datum<eT>::nan;
    }
  
@@ -116,10 +116,15 @@ op_imag::apply( Mat<typename T1::pod_type>& out, const mtOp<typename T1::pod_typ
  {
  arma_debug_sigprint();
  
-  typedef typename T1::pod_type T;
+  typedef typename T1::elem_type eT;
+  typedef typename T1::pod_type   T;
  
  const Proxy<T1> P(X.m);
  
+  if(is_cx<eT>::no)  { out.zeros(P.get_n_rows(), P.get_n_cols()); return; }
+  
+  // aliasing not possible at this point, as eT must be std::complex
+  
  const uword n_rows = P.get_n_rows();
  const uword n_cols = P.get_n_cols();
    
@@ -136,7 +141,7 @@ op_imag::apply( Mat<typename T1::pod_type>& out, const mtOp<typename T1::pod_typ
    
    for(uword i=0; i < n_elem; ++i)
      {
-      out_mem[i] = std::imag( A[i] );
+      out_mem[i] = access::tmp_imag( A[i] );
      }
    }
  else
@@ -144,7 +149,7 @@ op_imag::apply( Mat<typename T1::pod_type>& out, const mtOp<typename T1::pod_typ
    for(uword col=0; col < n_cols; ++col)
    for(uword row=0; row < n_rows; ++row)
      {
-      *out_mem = std::imag( P.at(row,col) );
+      *out_mem = access::tmp_imag( P.at(row,col) );
      out_mem++;
      }
    }
@@ -159,10 +164,15 @@ op_imag::apply( Cube<typename T1::pod_type>& out, const mtOpCube<typename T1::po
  {
  arma_debug_sigprint();
  
-  typedef typename T1::pod_type T;
+  typedef typename T1::elem_type eT;
+  typedef typename T1::pod_type   T;
  
  const ProxyCube<T1> P(X.m);
  
+  if(is_cx<eT>::no)  { out.zeros(P.get_n_rows(), P.get_n_cols(), P.get_n_slices()); return; }
+  
+  // aliasing not possible at this point, as eT must be std::complex
+  
  const uword n_rows   = P.get_n_rows();
  const uword n_cols   = P.get_n_cols();
  const uword n_slices = P.get_n_slices();
@@ -180,7 +190,7 @@ op_imag::apply( Cube<typename T1::pod_type>& out, const mtOpCube<typename T1::po
    
    for(uword i=0; i < n_elem; ++i)
      {
-      out_mem[i] = std::imag( A[i] );
+      out_mem[i] = access::tmp_imag( A[i] );
      }
    }
  else
@@ -189,7 +199,7 @@ op_imag::apply( Cube<typename T1::pod_type>& out, const mtOpCube<typename T1::po
    for(uword col=0;   col   < n_cols;   ++col  )
    for(uword row=0;   row   < n_rows;   ++row  )
      {
-      *out_mem = std::imag( P.at(row,col,slice) );
+      *out_mem = access::tmp_imag( P.at(row,col,slice) );
      out_mem++;
      }
    }
--- a/Show More
+++ b/Show More
Author	SHA1	Message	Date
conrad	5fe5a0303a	patch bump	2025-08-25 20:21:06 +10:00
conrad	58d012deaf	better wording	2025-07-12 18:46:57 +10:00
conrad	f459f8c05b	fix	2025-07-12 11:19:26 +10:00
conrad	f1e256cdc2	cleanup	2025-07-07 15:43:29 +10:00
conrad	9f9bc041a0	cleanup wordiing	2025-07-04 19:37:49 +10:00
conrad	5ab05fc6a1	fix	2025-06-30 12:42:43 +10:00
conrad	2c15e54506	fix	2025-06-12 12:07:00 +10:00
conrad	8ab2cf9589	patch bump	2025-05-21 15:14:53 +10:00
conrad	90f305d153	take into account possible inf and nan	2025-05-21 10:20:08 +10:00
conrad	4ddb01bfd7	restore previous attributes due to false positives in gcc 14.2	2025-05-20 13:46:33 +10:00
conrad	c4ffcac263	remove no longer needed typedef	2025-05-20 13:40:59 +10:00
conrad	83a38cc999	fix comment	2025-05-20 13:02:35 +10:00
conrad	5a80708b4a	optimisation for empty matrices	2025-05-20 12:46:37 +10:00
conrad	365af48910	take into account possible inf and nan	2025-05-20 12:33:58 +10:00
conrad	4a2cdf78e0	take into account possible inf and nan	2025-05-20 01:15:25 +10:00
conrad	e3f0ff95b9	minor layout fix	2025-05-19 12:49:25 +10:00
conrad	8cb64ddc3d	fix handling of corner-cases with nan and inf	2025-05-19 12:40:55 +10:00
conrad	fa1cc1925e	improve short summary for shift()	2025-05-18 16:24:52 +10:00
conrad	25f6093cef	conditionally undefine ARMA_IGNORE_DEPRECATED_MARKER	2025-05-15 13:33:25 +10:00
conrad	5d93b06328	better wording	2025-05-14 13:18:02 +10:00
conrad	210936459e	use .is_finite()	2025-05-13 16:20:20 +10:00
conrad	c36c50476a	add DOI	2025-05-12 23:31:03 +10:00
conrad	ef46b50c0c	deprecate stand-alone is_finite()	2025-05-07 12:56:51 +10:00
conrad	4d9a82294a	issue compile time warning for unused outputs	2025-05-02 14:14:07 +10:00
conrad	0ccf188ecd	use std:: prefix	2025-04-26 00:29:02 +10:00
conrad	caf55b8129	patch bump	2025-04-25 23:59:11 +10:00
conrad	70ad84b6fc	update pdf	2025-04-25 23:57:15 +10:00
conrad	413197f906	remove old papers	2025-04-25 23:55:18 +10:00
conrad	75a08f29d9	workaround for bugs in clang 20.1	2025-04-24 15:18:05 +10:00
conrad	10b1286069	refactor to use superlu::int_t	2025-04-23 16:27:05 +10:00
conrad	ff07666ca1	fix layout	2025-04-23 14:48:20 +10:00
conrad	f745a4fd64	embed icon image into html	2025-04-23 14:45:53 +10:00
conrad	e4dbea66f2	update paper links	2025-04-23 14:30:42 +10:00
conrad	c51a929a95	better wording	2025-04-23 14:28:26 +10:00
conrad	9384470685	better wording	2025-04-23 14:27:27 +10:00
conrad	df8c66b33d	bump min version of CMake	2025-04-23 14:26:51 +10:00
conrad	38db8d899e	add compile time warning	2025-04-23 14:08:11 +10:00
conrad	08bcbadd40	add compile time warning	2025-04-23 14:05:24 +10:00
conrad	f7d5691c0c	add support for superlu 7.0	2025-04-23 13:38:59 +10:00
conrad	c22443e423	rename	2025-04-23 13:33:39 +10:00
conrad	92a91a20d7	fix speed regression for inv of tiny matrices	2025-04-23 13:31:13 +10:00
conrad	c4a99aae01	add trace reduction	2025-04-23 13:25:35 +10:00
conrad	b613bdb306	fix calculation of s	2025-04-23 13:24:40 +10:00
conrad	f77f5c3d0f	simplifications	2025-04-11 11:50:11 +10:00
conrad	4e64c2e1b4	improve docs for element access via better ordering	2025-04-07 14:42:42 +10:00
conrad	2456f676ec	various spelling fixes	2025-03-25 23:17:48 +10:00
conrad	55560bbdf5	spelling fixes	2025-03-25 21:02:03 +10:00
conrad	6e880d8174	patch bump	2025-03-25 16:31:24 +10:00
conrad	0945718ee7	update pdf	2025-03-25 16:27:19 +10:00
conrad	cc20c82f81	rename pdf	2025-03-25 16:26:39 +10:00
conrad	12259d14c4	add mutex around fftw plan teardown	2025-03-20 22:54:43 +08:00
conrad	993fe423df	better wording	2025-03-19 23:52:25 +10:00
conrad	91970b89ff	better wording	2025-03-19 14:14:06 +10:00
conrad	5bcbfc2494	improve docs for vecnorm()	2025-03-19 14:03:51 +10:00
conrad	64204ba71e	minor simplifications	2025-03-18 19:58:20 +10:00
conrad	efbb9ed899	update citation details	2025-03-17 15:31:46 +10:00
conrad	79d8be32a4	update citation details	2025-03-10 15:31:46 +10:00
conrad	e5e59aad06	better wording	2025-03-10 12:49:05 +10:00
conrad	18836f6f54	minor fixes	2025-02-26 17:16:27 +10:00
conrad	12134220c6	add debug statements	2025-02-20 11:38:58 +10:00
conrad	54c8a6bce9	minor version bump	2025-02-17 14:23:00 +10:00
conrad	597779ad8e	remove commented out code	2025-02-17 14:03:24 +10:00
conrad	c7e0735820	mark as RC2	2025-02-15 15:20:49 +10:00
conrad	0a175fb070	comment out problematic workarounds	2025-02-15 11:56:07 +10:00
conrad	1ee1a6ea2d	remove problematic workaround	2025-02-15 11:49:54 +10:00
conrad	cb5686067d	minor speedup	2025-02-14 14:47:06 +10:00
conrad	5d92b91916	faster handling of sum(square(X))	2025-02-14 14:40:31 +10:00
conrad	dca21aa5f1	optimisation for spglue_schur	2025-02-14 14:32:11 +10:00
conrad	fce2e450b8	allow use of blas::dot()	2025-02-14 14:12:24 +10:00
conrad	28139a6da3	mark as RC1	2025-02-13 23:29:27 +10:00
conrad	5c45d78bba	workaroudns for old user code	2025-02-13 22:45:47 +10:00
conrad	0446d2eed8	minor layout fix	2025-02-13 22:45:12 +10:00
conrad	3f0e62abfa	workaround for old user code	2025-02-12 11:14:43 +10:00
conrad	da80bc561f	fix types	2025-02-12 10:04:02 +10:00
conrad	5a845df84b	fix type	2025-02-12 09:46:43 +10:00
conrad	c6219667df	link to local PDF	2025-02-12 09:46:18 +10:00
conrad	90b529091f	add arxiv pre-print	2025-02-11 13:08:36 +10:00
conrad	dc1ede314e	update citation	2025-02-10 15:34:37 +10:00
conrad	fcd466d213	update copyright year	2025-02-10 15:29:46 +10:00
conrad	8ac6a694f0	update list of changes	2025-02-10 15:00:36 +10:00
conrad	9e6bda097d	pre-calculate abs(cx) for faster sorting of complex matrices	2025-02-10 14:44:19 +10:00
conrad	0bea1b2f4d	simplify comment	2025-02-10 12:03:26 +10:00
conrad	08a15ab6d2	simplifications	2025-02-07 15:44:24 +10:00
conrad	8f76c2f5a7	use temporary matrix to explicitly avoid aliasing	2025-02-07 15:32:51 +10:00
conrad	eea1835183	specialised handling for pow(x,0.5)	2025-02-07 12:07:38 +10:00
conrad	8e02066f52	better wording	2025-02-07 11:51:28 +10:00
conrad	cc5d113370	reinterpret pow(x,0.5) as sqrt(x)	2025-02-07 11:32:05 +10:00
conrad	8903fd84a5	restore selective use_mp for pow()	2025-02-07 11:15:52 +10:00
conrad	08d3312192	reinterpret pow(x,0.5) as sqrt(x)	2025-02-07 11:06:51 +10:00
conrad	9944415e12	minor fix	2025-02-06 14:44:15 +10:00
conrad	beb4ee0e39	remove unused variable	2025-02-06 14:44:00 +10:00
conrad	4ab6a53ab5	minor layout fix	2025-02-06 14:40:31 +10:00
conrad	99b11428df	use std:: prefix	2025-02-05 16:34:40 +10:00
conrad	955abfeeab	refactor to avoid using union	2025-02-05 16:24:14 +10:00
conrad	e5b1d0bd70	replace weird typedefs with C++11 type aliases	2025-02-05 13:01:33 +10:00
conrad	b70e88d8ad	minor fix	2025-02-04 13:40:24 +10:00
conrad	e7ecee425f	remove checking for matrix mul chains	2025-02-03 14:55:59 +10:00
conrad	7004caa4b4	handle more corner cases	2025-02-03 14:41:16 +10:00
conrad	22bbf2ad79	better debug message	2025-02-03 14:10:20 +10:00
conrad	73294812c3	rename unwrap variables	2025-02-03 13:22:42 +10:00
conrad	52b372f7f0	remove partial_unwrap_check	2025-02-03 13:15:07 +10:00
conrad	1f61355864	refactor to use partial_unwrap	2025-02-03 13:11:24 +10:00
conrad	95ac9ebde7	refactor to use partial_unwrap	2025-02-03 12:56:34 +10:00
conrad	b244199997	print internal flags	2025-02-03 11:23:40 +10:00
conrad	c689776db2	preserve more type information	2025-02-02 23:52:17 +10:00
conrad	37ca4fe05c	update list of changes	2025-02-02 11:22:44 +10:00
conrad	4a299aa8ec	update list of changes	2025-01-28 14:28:05 +10:00
conrad	e0162f0a84	take into account use_mp flag	2025-01-28 14:24:53 +10:00
conrad	cb9c78a9e3	faster handling of square() and pow(x,2)	2025-01-28 12:25:17 +10:00
conrad	de479975af	simplifications	2025-01-28 11:05:11 +10:00
conrad	d57ad2a29a	optimisation for spop_square	2025-01-27 22:47:34 +10:00
conrad	42a35f8cc3	optimisation for eop_square	2025-01-27 20:36:51 +10:00
conrad	a75737fef7	use quasi_unwrap earlier	2025-01-27 00:10:13 +10:00
conrad	426ce1372a	directly avoid Proxy if possible	2025-01-25 22:29:44 +10:00
conrad	bce19ebcbd	directly avoid Proxy if possible	2025-01-25 22:00:08 +10:00
conrad	59477c0dbc	restore previous version	2025-01-24 13:08:05 +10:00
conrad	bec851861b	harmonise to use proxy_type	2025-01-24 11:43:39 +10:00
conrad	4e2f157604	reinterpret pow(x,2) as square(x)	2025-01-24 10:56:07 +10:00
conrad	2b513d3136	add TODO	2025-01-23 18:09:05 +10:00
conrad	efac5b3148	reinterpret pow(x,2) as square(x)	2025-01-23 17:11:42 +10:00
conrad	daa5deedf3	explicitly handle another corner case	2025-01-22 23:50:27 +10:00
conrad	e658f0db17	restore old option for now	2025-01-22 23:23:02 +10:00
conrad	bcee6924b7	update list of changes	2025-01-20 14:14:37 +10:00
conrad	089d2b965d	unwrap if openmp can be used	2025-01-20 00:15:23 +10:00
conrad	ee9f42faff	simplifications	2025-01-19 23:41:41 +10:00
conrad	92f4bd5a24	simplifications	2025-01-19 23:01:35 +10:00
conrad	25383133c8	avoid using the proxy if we have direct memory access	2025-01-19 22:47:16 +10:00
conrad	bde6942425	update list of changes	2025-01-19 00:17:26 +10:00
conrad	b9461be144	faster handling of single submatirx columns	2025-01-18 23:58:31 +10:00
conrad	c83caaa0da	more efficient handling of expressions	2025-01-18 23:56:57 +10:00
conrad	0ec6b425e6	avoid using the proxy if we have direct memory access	2025-01-18 23:16:11 +10:00
conrad	785bfcb086	avoid using the proxy if we have direct memory access	2025-01-18 23:03:49 +10:00
conrad	bb515b8a1f	simplifications	2025-01-17 21:33:16 +10:00
conrad	cb57dc3815	generate zeros quickly for non-complex inputs	2025-01-17 21:24:45 +10:00
conrad	f474274d92	preserve type information	2025-01-17 21:19:04 +10:00
conrad	1629b41828	add warning	2025-01-17 14:08:13 +10:00
conrad	20afeff672	fix	2025-01-17 14:07:55 +10:00
conrad	78694113e3	cleanup	2025-01-17 13:41:00 +10:00
conrad	acd2a5fcad	optimisations for common cases	2025-01-17 13:37:58 +10:00
conrad	32fd3716fe	optimisations for common cases	2025-01-17 13:33:13 +10:00
conrad	64de62450d	exploit known dimensions	2025-01-17 13:13:12 +10:00
conrad	81e7d78488	don't fill with zeros as the memory will be immediately overwritten	2025-01-17 13:05:54 +10:00
conrad	b8c402b663	avoid copying memory if possible	2025-01-17 13:01:07 +10:00
conrad	328f0fe573	exploit known dimensions	2025-01-17 12:51:38 +10:00
conrad	27db1972be	fix for cornercase interaction with .slice()	2025-01-17 12:32:21 +10:00
conrad	b0fe0d3cff	simplifications	2025-01-17 12:23:08 +10:00
conrad	522cfdce58	fix	2025-01-16 16:34:06 +10:00
conrad	dbf64f011c	explicitly check for allowed transformations instead of implicit via n_elem	2025-01-16 16:29:13 +10:00
conrad	f90249082c	simplification	2025-01-16 16:14:14 +10:00
conrad	3b5fb76081	update list of changes	2025-01-16 16:04:47 +10:00
conrad	af912a6177	fix	2025-01-16 15:59:42 +10:00
conrad	803c589b17	expand field class with initial implementations of .reshape() and .resize()	2025-01-16 15:50:59 +10:00
conrad	5eb6b5b5fc	fix comments	2025-01-12 23:13:12 +10:00
conrad	db71a4affb	simplifications	2025-01-12 22:57:53 +10:00
conrad	0936110fbc	update list of changes	2025-01-09 14:54:02 +10:00
conrad	f905ee4b4a	reinterpret pow(x,2) as square(x)	2025-01-09 14:39:11 +10:00
conrad	7ff62d22b0	simplifications	2025-01-09 14:30:55 +10:00
conrad	c004318b0f	fix	2025-01-09 14:25:20 +10:00
conrad	3632831126	reinterpret pow(x,2) as square(x)	2025-01-09 14:24:14 +10:00
conrad	ebc6f94fcb	simplification	2025-01-09 14:18:42 +10:00
conrad	b4ede23c5c	reinterpret pow(x,2) as square(x)	2025-01-09 14:17:15 +10:00
conrad	feca3eae7f	fix	2025-01-09 12:50:41 +10:00
conrad	0a76ab431b	remove support for ARMA_IGNORE_DEPRECATED_MARKER	2025-01-02 15:09:00 +10:00
conrad	24b059344e	new section for list of changes	2025-01-02 13:11:14 +10:00
conrad	ceaa020a33	better layout for sinc()	2025-01-02 13:07:25 +10:00
conrad	f9e308b2fd	improve docs for toeplitz()	2025-01-02 00:02:21 +10:00
conrad	153c1363ea	simplify docs for toeplitz()	2025-01-01 23:28:01 +10:00
conrad	316c7e4abf	simplification; alias checks not required	2024-12-25 00:22:08 +10:00
conrad	469cb91a74	reduce stack usage by reducing unnecessary alignment markers	2024-12-18 14:08:42 +10:00
conrad	77cc3e9287	fix prototypes for xSYCON and xHECON	2024-12-17 12:17:02 +10:00
conrad	e037bd9a5e	reduce bloat	2024-12-16 13:58:00 +10:00
conrad	7138865c14	add note	2024-12-10 12:05:52 +10:00
conrad	8f94af0afb	explicitly qualify inherited variables and functions	2024-12-04 12:24:44 +10:00
conrad	dfdbae2570	explicitly qualify inherited variables and functions	2024-12-04 12:16:50 +10:00
conrad	11e69548f2	preserve vector type information	2024-12-03 16:25:38 +10:00
conrad	6f6cc83862	fix	2024-12-03 16:19:52 +10:00
conrad	0d2c5de574	preserve vector type information	2024-12-02 22:44:30 +10:00
conrad	07e5b6d1d6	avoid clashes with C++20 std::span	2024-12-02 16:23:54 +10:00
conrad	df8741c86b	explicitly qualify inherited variables and functions	2024-12-02 16:13:48 +10:00
conrad	ffb9ef3c52	more selective detection of symmetric/hermitian matrices	2024-11-28 16:15:56 +11:00
conrad	b87cc838a3	more selective detection of symmetric/hermitian matrices	2024-11-27 11:56:43 +11:00
conrad	9273614f58	fix	2024-11-24 17:12:41 +10:00
conrad	cdedada64a	revert	2024-11-23 16:40:29 +10:00
conrad	b7cf2fd503	remove functions from Mat that are present in Base	2024-11-23 13:49:04 +10:00
conrad	b3329329a7	expand deprecation warnings with messages	2024-11-23 00:42:36 +10:00
conrad	20bb9d8bcb	ensure the first index of extremum value is found	2024-11-22 12:17:20 +10:00
conrad	6ecbed3813	remove unused code	2024-11-20 17:41:01 +10:00
conrad	123dd0108d	use partial_unwrap instead of partial_unwrap_check	2024-11-20 17:26:49 +10:00
conrad	79793a7a97	mark as unstable	2024-11-20 12:27:28 +10:00
conrad	a2c8b74a35	add TODOs	2024-11-20 11:45:56 +10:00