version bump

README.md

@@ -38,8 +38,7 @@ Copyright 2017-2025 Data61 / CSIRO
 13. [API Stability and Version Policy](#13-api-stability-and-version-policy)
 14. [Bug Reports and Frequently Asked Questions](#14-bug-reports-and-frequently-asked-questions)
 
-15. [MEX Interface to Octave/Matlab](#15-mex-interface-to-octavematlab)
-16. [Related Software Using Armadillo](#16-related-software-using-armadillo)
+15. [Related Software Using Armadillo](#15-related-software-using-armadillo)
 
 ---
 
@@ -452,14 +451,7 @@ https://arma.sourceforge.net/faq.html
 
 ---
 
-### 15: MEX Interface to Octave/Matlab
-
-The `mex_interface` folder contains examples of how to interface
-Octave/Matlab with C++ code that uses Armadillo matrices.
-
----
-
-### 16: Related Software Using Armadillo
+### 15: Related Software Using Armadillo
 
 * MLPACK: extensive library of machine learning algorithms  
   https://mlpack.org

docs.html

@@ -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 15.0</b></big>
+<big><b>API Documentation for Armadillo 15.2</b></big>
 <br>
 <br>
 <br>
@@ -360,10 +360,11 @@ Conrad Sanderson and Ryan Curtin.
 <tr><td><a href="#zeros_standalone">zeros</a></td><td>&nbsp;</td><td>generate object filled with zeros</td></tr>
 <tr style="background-color: #F5F5F5;"><td><a href="#randu">randu</a></td><td>&nbsp;</td><td>generate object with random values (uniform distribution)</td></tr>
 <tr style="background-color: #F5F5F5;"><td><a href="#randn">randn</a></td><td>&nbsp;</td><td>generate object with random values (normal distribution)</td></tr>
-<tr style="background-color: #F5F5F5;"><td><a href="#randg">randg</a></td><td>&nbsp;</td><td>generate object with random values (gamma distribution)</td></tr>
+<tr style="background-color: #F5F5F5;"><td><a href="#rande">rande</a></td><td>&nbsp;</td><td>generate object with random values (exponential distribution)</td></tr>
+<tr><td><a href="#randg">randg</a></td><td>&nbsp;</td><td>generate object with random values (gamma distribution)</td></tr>
 <tr><td><a href="#randi">randi</a></td><td>&nbsp;</td><td>generate object with random integer values in specified interval</td></tr>
 <tr><td><a href="#speye">speye</a></td><td>&nbsp;</td><td>generate sparse identity matrix</td></tr>
-<tr><td><a href="#spones">spones</a></td><td>&nbsp;</td><td>generate sparse matrix with non-zero elements set to one</td></tr>
+<tr style="background-color: #F5F5F5;"><td><a href="#spones">spones</a></td><td>&nbsp;</td><td>generate sparse matrix with non-zero elements set to one</td></tr>
 <tr style="background-color: #F5F5F5;"><td><a href="#sprandu_sprandn">sprandu&nbsp;/&nbsp;sprandn</a></td><td>&nbsp;</td><td>generate sparse matrix with non-zero elements set to random values</td></tr>
 <tr style="background-color: #F5F5F5;"><td><a href="#toeplitz">toeplitz</a></td><td>&nbsp;</td><td>generate Toeplitz matrix</td></tr>
 </tbody>
@@ -384,59 +385,59 @@ Conrad Sanderson and Ryan Curtin.
 <tr><td><a href="#arg">arg</a></td><td>&nbsp;</td><td>phase angle of each element</td></tr>
 <tr style="background-color: #F5F5F5;"><td><a href="#as_scalar">as_scalar</a></td><td>&nbsp;</td><td>convert 1x1 matrix to pure scalar</td></tr>
 <tr style="background-color: #F5F5F5;"><td><a href="#balance">balance</a></td><td>&nbsp;</td><td>balance matrix</td></tr>
-<tr style="background-color: #F5F5F5;"><td><a href="#clamp">clamp</a></td><td>&nbsp;</td><td>obtain clamped elements according to given limits</td></tr>
+<tr style="background-color: #F5F5F5;"><td><a href="#circshift">circshift</a></td><td>&nbsp;</td><td>circular shift of elements</td></tr>
+<tr><td><a href="#clamp">clamp</a></td><td>&nbsp;</td><td>obtain clamped elements according to given limits</td></tr>
 <tr><td><a href="#cond">cond</a></td><td>&nbsp;</td><td>condition number of matrix</td></tr>
 <tr><td><a href="#conj">conj</a></td><td>&nbsp;</td><td>obtain complex conjugate of each element</td></tr>
-<tr><td><a href="#conv_to">conv_to</a></td><td>&nbsp;</td><td>convert/cast between matrix types</td></tr>
+<tr style="background-color: #F5F5F5;"><td><a href="#conv_to">conv_to</a></td><td>&nbsp;</td><td>convert/cast between matrix types</td></tr>
 <tr style="background-color: #F5F5F5;"><td><a href="#cross">cross</a></td><td>&nbsp;</td><td>cross product</td></tr>
 <tr style="background-color: #F5F5F5;"><td><a href="#cumsum">cumsum</a></td><td>&nbsp;</td><td>cumulative sum</td></tr>
-<tr style="background-color: #F5F5F5;"><td><a href="#cumprod">cumprod</a></td><td>&nbsp;</td><td>cumulative product</td></tr>
+<tr><td><a href="#cumprod">cumprod</a></td><td>&nbsp;</td><td>cumulative product</td></tr>
 <tr><td><a href="#det">det</a></td><td>&nbsp;</td><td>determinant</td></tr>
 <tr><td><a href="#diagmat">diagmat</a></td><td>&nbsp;</td><td>generate diagonal matrix from given matrix or vector</td></tr>
-<tr><td><a href="#diagvec">diagvec</a></td><td>&nbsp;</td><td>extract specified diagonal</td></tr>
+<tr style="background-color: #F5F5F5;"><td><a href="#diagvec">diagvec</a></td><td>&nbsp;</td><td>extract specified diagonal</td></tr>
 <tr style="background-color: #F5F5F5;"><td><a href="#diags_spdiags">diags&nbsp;/&nbsp;spdiags</a></td><td>&nbsp;</td><td>generate band matrix from given set of vectors</td></tr>
 <tr style="background-color: #F5F5F5;"><td><a href="#diff">diff</a></td><td>&nbsp;</td><td>differences between adjacent elements</td></tr>
-<tr style="background-color: #F5F5F5;"><td><a href="#dot">dot&nbsp;/&nbsp;cdot&nbsp;/&nbsp;norm_dot</a></td><td>&nbsp;</td><td>dot product</td></tr>
+<tr><td><a href="#dot">dot&nbsp;/&nbsp;cdot&nbsp;/&nbsp;norm_dot</a></td><td>&nbsp;</td><td>dot product</td></tr>
 <tr><td><a href="#eps">eps</a></td><td>&nbsp;</td><td>obtain distance of each element to next largest floating point representation</td></tr>
 <tr><td><a href="#expmat">expmat</a></td><td>&nbsp;</td><td>matrix exponential</td></tr>
-<tr><td><a href="#expmat_sym">expmat_sym</a></td><td>&nbsp;</td><td>matrix exponential of symmetric matrix</td></tr>
+<tr style="background-color: #F5F5F5;"><td><a href="#expmat_sym">expmat_sym</a></td><td>&nbsp;</td><td>matrix exponential of symmetric matrix</td></tr>
 <tr style="background-color: #F5F5F5;"><td><a href="#find">find</a></td><td>&nbsp;</td><td>find indices of non-zero elements, or elements satisfying a relational condition</td></tr>
 <tr style="background-color: #F5F5F5;"><td><a href="#find_finite_nonfinite">find_finite / nonfinite</a></td><td>&nbsp;</td><td>find indices of finite / non-finite elements</td></tr>
-<tr style="background-color: #F5F5F5;"><td><a href="#find_nan_nonnan">find_nan / nonnan</a></td><td>&nbsp;</td><td>find indices of NaN / non-NaN elements</td></tr>
+<tr><td><a href="#find_nan_nonnan">find_nan / nonnan</a></td><td>&nbsp;</td><td>find indices of NaN / non-NaN elements</td></tr>
 <tr><td><a href="#find_unique">find_unique</a></td><td>&nbsp;</td><td>find indices of unique elements</td></tr>
 <tr><td><a href="#flip">fliplr&nbsp;/&nbsp;flipud</a></td><td>&nbsp;</td><td>flip matrix left to right or upside down</td></tr>
-<tr><td><a href="#imag_real">imag / real</a></td><td>&nbsp;</td><td>extract imaginary/real part</td></tr>
+<tr style="background-color: #F5F5F5;"><td><a href="#imag_real">imag / real</a></td><td>&nbsp;</td><td>extract imaginary/real part</td></tr>
 <tr style="background-color: #F5F5F5;"><td><a href="#ind2sub">ind2sub</a></td><td>&nbsp;</td><td>convert linear index to subscripts</td></tr>
 <tr style="background-color: #F5F5F5;"><td><a href="#index_min_and_index_max_standalone">index_min / index_max</a></td><td>&nbsp;</td><td>indices of extremum values</td></tr>
-<tr style="background-color: #F5F5F5;"><td><a href="#inplace_trans">inplace_trans</a></td><td>&nbsp;</td><td>in-place transpose</td></tr>
+<tr><td><a href="#inplace_trans">inplace_trans</a></td><td>&nbsp;</td><td>in-place transpose</td></tr>
 <tr><td><a href="#intersect">intersect</a></td><td>&nbsp;</td><td>find common elements in two vectors/matrices</td></tr>
 <tr><td><a href="#join">join_rows / join_cols</a></td><td>&nbsp;</td><td>concatenation of matrices</td></tr>
-<tr><td><a href="#join_slices">join_slices</a></td><td>&nbsp;</td><td>concatenation of cubes</td></tr>
+<tr style="background-color: #F5F5F5;"><td><a href="#join_slices">join_slices</a></td><td>&nbsp;</td><td>concatenation of cubes</td></tr>
 <tr style="background-color: #F5F5F5;"><td><a href="#kron">kron</a></td><td>&nbsp;</td><td>Kronecker tensor product</td></tr>
 <tr style="background-color: #F5F5F5;"><td><a href="#log_det">log_det</a></td><td>&nbsp;</td><td>log determinant</td></tr>
-<tr style="background-color: #F5F5F5;"><td><a href="#log_det_sympd">log_det_sympd</a></td><td>&nbsp;</td><td>log determinant of symmetric positive definite matrix</td></tr>
+<tr><td><a href="#log_det_sympd">log_det_sympd</a></td><td>&nbsp;</td><td>log determinant of symmetric positive definite matrix</td></tr>
 <tr><td><a href="#logmat">logmat</a></td><td>&nbsp;</td><td>matrix logarithm</td></tr>
 <tr><td><a href="#logmat_sympd">logmat_sympd</a></td><td>&nbsp;</td><td>matrix logarithm of symmetric matrix</td></tr>
-<tr><td><a href="#min_and_max">min&nbsp;/&nbsp;max</a></td><td>&nbsp;</td><td>return extremum values</td></tr>
+<tr style="background-color: #F5F5F5;"><td><a href="#min_and_max">min&nbsp;/&nbsp;max</a></td><td>&nbsp;</td><td>return extremum values</td></tr>
 <tr style="background-color: #F5F5F5;"><td><a href="#norm">norm</a></td><td>&nbsp;</td><td>various norms of vectors and matrices</td></tr>
 <tr style="background-color: #F5F5F5;"><td><a href="#norm2est">norm2est</a></td><td>&nbsp;</td><td>fast estimate of the matrix 2-norm</td></tr>
-<tr style="background-color: #F5F5F5;"><td><a href="#normalise">normalise</a></td><td>&nbsp;</td><td>normalise vectors to unit <i>p</i>-norm</td></tr>
+<tr><td><a href="#normalise">normalise</a></td><td>&nbsp;</td><td>normalise vectors to unit <i>p</i>-norm</td></tr>
 <tr><td><a href="#nonzeros">nonzeros</a></td><td>&nbsp;</td><td>extract all non-zero elements</td></tr>
 <tr><td><a href="#omit_nan_nonfinite">omit_nan / nonfinite</a></td><td>&nbsp;</td><td>extract all elements that are non-NaN / only finite</td></tr>
-<tr><td><a href="#pow">pow</a></td><td>&nbsp;</td><td>element-wise power</td></tr>
+<tr style="background-color: #F5F5F5;"><td><a href="#pow">pow</a></td><td>&nbsp;</td><td>element-wise power</td></tr>
 <tr style="background-color: #F5F5F5;"><td><a href="#powmat">powmat</a></td><td>&nbsp;</td><td>matrix power</td></tr>
 <tr style="background-color: #F5F5F5;"><td><a href="#prod">prod</a></td><td>&nbsp;</td><td>product of elements</td></tr>
-<tr style="background-color: #F5F5F5;"><td><a href="#rank">rank</a></td><td>&nbsp;</td><td>rank of matrix</td></tr>
+<tr><td><a href="#rank">rank</a></td><td>&nbsp;</td><td>rank of matrix</td></tr>
 <tr><td><a href="#rcond">rcond</a></td><td>&nbsp;</td><td>reciprocal condition number</td></tr>
 <tr><td><a href="#repelem">repelem</a></td><td>&nbsp;</td><td>replicate elements</td></tr>
-<tr><td><a href="#replace_standalone">replace</a></td><td>&nbsp;</td><td>replace specific elements with a new value</td></tr>
+<tr style="background-color: #F5F5F5;"><td><a href="#replace_standalone">replace</a></td><td>&nbsp;</td><td>replace specific elements with a new value</td></tr>
 <tr style="background-color: #F5F5F5;"><td><a href="#repmat">repmat</a></td><td>&nbsp;</td><td>replicate matrix in block-like fashion</td></tr>
 <tr style="background-color: #F5F5F5;"><td><a href="#repcube">repcube</a></td><td>&nbsp;</td><td>replicate cube in block-like fashion</td></tr>
-<tr style="background-color: #F5F5F5;"><td><a href="#reshape">reshape</a></td><td>&nbsp;</td><td>change size while keeping elements</td></tr>
+<tr><td><a href="#reshape">reshape</a></td><td>&nbsp;</td><td>change size while keeping elements</td></tr>
 <tr><td><a href="#resize">resize</a></td><td>&nbsp;</td><td>change size while keeping elements and preserving layout</td></tr>
 <tr><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 style="background-color: #F5F5F5;"><td><a href="#shift">shift</a></td><td>&nbsp;</td><td>circular shift of elements</td></tr>
+<tr style="background-color: #F5F5F5;"><td><a href="#roots">roots</a></td><td>&nbsp;</td><td>roots of polynomial</td></tr>
 <tr style="background-color: #F5F5F5;"><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><td><a href="#sort">sort</a></td><td>&nbsp;</td><td>sort elements</td></tr>
@@ -701,16 +702,6 @@ For convenience, the following matrix typedefs are defined:
       <code>Mat&lt;<a href="#uword">sword</a>&gt;</code>
       </td>
     </tr>
-  </tbody>
-</table>
-</li>
-<br>
-<li>
-If supported by compiler and hardware, additional typedefs for matrices with half-precision element types are defined as: 
-<br>
-<br>
-<table style="text-align: left;" border="0" cellpadding="2" cellspacing="2">
-  <tbody>
     <tr>
       <td style="vertical-align: top; text-align: right;">
       <code>hmat</code>
@@ -721,6 +712,9 @@ If supported by compiler and hardware, additional typedefs for matrices with hal
       <td style="vertical-align: top;">
       <code>Mat&lt;<a href="#fp16_type">fp16</a>&gt;</code>
       </td>
+      <td style="vertical-align: top;">
+      &nbsp;(only available if supported by compiler and hardware)
+      </td>
     </tr>
     <tr>
       <td style="vertical-align: top; text-align: right;">
@@ -732,6 +726,9 @@ If supported by compiler and hardware, additional typedefs for matrices with hal
       <td style="vertical-align: top;">
       <code>Mat&lt;<a href="#fp16_type">cx_fp16</a>&gt;</code>
       </td>
+      <td style="vertical-align: top;">
+      &nbsp;(only available if supported by compiler and hardware)
+      </td>
     </tr>
   </tbody>
 </table>
@@ -776,7 +773,7 @@ Constructors:
 </li>
 <br>
 <li>
-The elements can be explicitly initialised during construction by specifying <i>fill_form</i>,
+Matrix elements can be explicitly initialised during construction by specifying <i>fill_form</i>,
 which is one of:
 <ul>
 <table>
@@ -1095,20 +1092,6 @@ For convenience, the following column vector typedefs are defined:
       <code>Col&lt;<a href="#uword">sword</a>&gt;</code>
       </td>
     </tr>
-  </tbody>
-</table>
-</li>
-<br>
-<li>
-The <b><i>vec</i></b> and <b><i>colvec</i></b> types have the <b>same meaning</b> and are used <b>interchangeably</b>
-</li>
-<br>
-<li>
-If supported by compiler and hardware, additional typedefs for column vectors with half-precision element types are defined as: 
-<br>
-<br>
-<table style="text-align: left;" border="0" cellpadding="2" cellspacing="2">
-  <tbody>
     <tr>
       <td style="vertical-align: top; text-align: right;">
       <code>hvec</code>
@@ -1125,6 +1108,9 @@ If supported by compiler and hardware, additional typedefs for column vectors wi
       <td style="vertical-align: top;">
       <code>Col&lt;<a href="#fp16_type">fp16</a>&gt;</code>
       </td>
+      <td style="vertical-align: top;">
+      &nbsp;(only available if supported by compiler and hardware)
+      </td>
     </tr>
     <tr>
       <td style="vertical-align: top; text-align: right;">
@@ -1142,12 +1128,19 @@ If supported by compiler and hardware, additional typedefs for column vectors wi
       <td style="vertical-align: top;">
       <code>Col&lt;<a href="#fp16_type">cx_fp16</a>&gt;</code>
       </td>
+      <td style="vertical-align: top;">
+      &nbsp;(only available if supported by compiler and hardware)
+      </td>
     </tr>
   </tbody>
 </table>
 </li>
 <br>
 <li>
+The <b><i>vec</i></b> and <b><i>colvec</i></b> types have the <b>same meaning</b> and are used <b>interchangeably</b>
+</li>
+<br>
+<li>
 In this documentation, the types <i>vec</i> or <i>colvec</i> are used for convenience; it is possible to use other column vector types instead, eg. <i>fvec</i>, <i>fcolvec</i>
 </li>
 <br>
@@ -1178,11 +1171,7 @@ Constructors:
 </li>
 <br>
 <li>
-<b>Caveat:</b>
-<ul>
-<li>in Armadillo &gt;= 10.5, the default initialisation is <code>fill::zeros</code></li>
-<li>in Armadillo &lt;= 10.4, the default initialisation is <code>fill::none</code></li>
-</ul>
+Vector elements can be explicitly initialised during construction by specifying <i>fill_form</i>, as per the <a href="#Mat">Mat</a> class
 </li>
 <br>
 <a name="adv_constructors_col"></a>
@@ -1390,16 +1379,6 @@ For convenience, the following row vector typedefs are defined:
       <code>Row&lt;<a href="#uword">sword</a>&gt;</code>
       </td>
     </tr>
-  </tbody>
-</table>
-</li>
-<br>
-<li>
-If supported by compiler and hardware, additional typedefs for row vectors with half-precision element types are defined as: 
-<br>
-<br>
-<table style="text-align: left;" border="0" cellpadding="2" cellspacing="2">
-  <tbody>
     <tr>
       <td style="vertical-align: top; text-align: right;">
       <code>hrowvec</code>
@@ -1410,6 +1389,9 @@ If supported by compiler and hardware, additional typedefs for row vectors with
       <td style="vertical-align: top;">
       <code>Row&lt;<a href="#fp16_type">fp16</a>&gt;</code>
       </td>
+      <td style="vertical-align: top;">
+      &nbsp;(only available if supported by compiler and hardware)
+      </td>
     </tr>
     <tr>
       <td style="vertical-align: top; text-align: right;">
@@ -1421,6 +1403,9 @@ If supported by compiler and hardware, additional typedefs for row vectors with
       <td style="vertical-align: top;">
       <code>Row&lt;<a href="#fp16_type">cx_fp16</a>&gt;</code>
       </td>
+      <td style="vertical-align: top;">
+      &nbsp;(only available if supported by compiler and hardware)
+      </td>
     </tr>
   </tbody>
 </table>
@@ -1458,11 +1443,7 @@ Constructors:
 </li>
 <br>
 <li>
-<b>Caveat:</b>
-<ul>
-<li>in Armadillo &gt;= 10.5, the default initialisation is <code>fill::zeros</code></li>
-<li>in Armadillo &lt;= 10.4, the default initialisation is <code>fill::none</code></li>
-</ul>
+Vector elements can be explicitly initialised during construction by specifying <i>fill_form</i>, as per the <a href="#Mat">Mat</a> class
 </li>
 <br>
 <a name="adv_constructors_row"></a>
@@ -1680,16 +1661,6 @@ For convenience, the following cube typedefs are defined:
       <code>Cube&lt;<a href="#uword">sword</a>&gt;</code>
       </td>
     </tr>
-  </tbody>
-</table>
-</li>
-<br>
-<li>
-If supported by compiler and hardware, additional typedefs for cubes with half-precision element types are defined as: 
-<br>
-<br>
-<table style="text-align: left;" border="0" cellpadding="2" cellspacing="2">
-  <tbody>
     <tr>
       <td style="vertical-align: top; text-align: right;">
       <code>hcube</code>
@@ -1700,6 +1671,9 @@ If supported by compiler and hardware, additional typedefs for cubes with half-p
       <td style="vertical-align: top;">
       <code>Cube&lt;<a href="#fp16_type">fp16</a>&gt;</code>
       </td>
+      <td style="vertical-align: top;">
+      &nbsp;(only available if supported by compiler and hardware)
+      </td>
     </tr>
     <tr>
       <td style="vertical-align: top; text-align: right;">
@@ -1711,6 +1685,9 @@ If supported by compiler and hardware, additional typedefs for cubes with half-p
       <td style="vertical-align: top;">
       <code>Cube&lt;<a href="#fp16_type">cx_fp16</a>&gt;</code>
       </td>
+      <td style="vertical-align: top;">
+      &nbsp;(only available if supported by compiler and hardware)
+      </td>
     </tr>
   </tbody>
 </table>
@@ -1740,7 +1717,7 @@ Constructors:
 </li>
 <br>
 <li>
-The elements can be explicitly initialised during construction by specifying <i>fill_form</i>,
+Cube elements can be explicitly initialised during construction by specifying <i>fill_form</i>,
 which is one of:
 <ul>
 <table>
@@ -1966,7 +1943,8 @@ The root sparse matrix class is <b>SpMat&lt;</b><i>type</i><b>&gt;</b>, where <i
 <br>
 <li>
 For convenience, the following sparse matrix typedefs are defined:
-<ul>
+<br>
+<br>
 <table style="text-align: left;" border="0" cellpadding="2" cellspacing="2">
   <tbody>
     <tr>
@@ -2057,17 +2035,6 @@ For convenience, the following sparse matrix typedefs are defined:
       <code>SpMat&lt;<a href="#uword">sword</a>&gt;</code>
       </td>
     </tr>
-  </tbody>
-</table>
-</ul>
-</li>
-<br>
-<li>
-If supported by compiler and hardware, additional typedefs for sparse matrices with half-precision element types are defined as: 
-<br>
-<br>
-<table style="text-align: left;" border="0" cellpadding="2" cellspacing="2">
-  <tbody>
     <tr>
       <td style="vertical-align: top; text-align: right;">
       <code>sp_hmat</code>
@@ -2078,6 +2045,9 @@ If supported by compiler and hardware, additional typedefs for sparse matrices w
       <td style="vertical-align: top;">
       <code>SpMat&lt;<a href="#fp16_type">fp16</a>&gt;</code>
       </td>
+      <td style="vertical-align: top;">
+      &nbsp;(only available if supported by compiler and hardware)
+      </td>
     </tr>
     <tr>
       <td style="vertical-align: top; text-align: right;">
@@ -2089,6 +2059,9 @@ If supported by compiler and hardware, additional typedefs for sparse matrices w
       <td style="vertical-align: top;">
       <code>SpMat&lt;<a href="#fp16_type">cx_fp16</a>&gt;</code>
       </td>
+      <td style="vertical-align: top;">
+      &nbsp;(only available if supported by compiler and hardware)
+      </td>
     </tr>
   </tbody>
 </table>
@@ -2200,7 +2173,7 @@ The following subset of operations &amp; functions is available for sparse matri
 <li>element-wise functions: <a href="#abs">abs()</a>, <a href="#misc_fns">cbrt()</a>, <a href="#misc_fns">ceil()</a>, <a href="#conj">conj()</a>, <a href="#misc_fns">floor()</a>, <a href="#imag_real">imag()</a>, <a href="#imag_real">real()</a>, <a href="#misc_fns">round()</a>, <a href="#misc_fns">sign()</a>, <a href="#misc_fns">sqrt()</a>, <a href="#misc_fns">square()</a>, <a href="#misc_fns">trunc()</a></li>
 <li>scalar functions of matrices: <a href="#accu">accu()</a>, <a href="#as_scalar">as_scalar()</a>, <a href="#dot">dot()</a>, <a href="#norm">norm()</a>, <a href="#norm2est">norm2est()</a>, <a href="#trace">trace()</a></li>
 <li>vector valued functions of matrices: <a href="#diagvec">diagvec()</a>, <a href="#min_and_max">min()</a>, <a href="#min_and_max">max()</a>, <a href="#nonzeros">nonzeros()</a>, <a href="#sum">sum()</a>, <a href="#stats_fns">mean()</a>, <a href="#stats_fns">stddev()</a>, <a href="#stats_fns">var()</a>, <a href="#vecnorm">vecnorm()</a>, <a href="#vectorise">vectorise()</a></li>
-<li>matrix valued functions of matrices: <a href="#clamp">clamp()</a>, <a href="#diagmat">diagmat()</a>, <a href="#diags_spdiags">spdiags()</a>, <a href="#flip">flipud()/fliplr()</a>, <a href="#join">join_rows()</a>, <a href="#join">join_cols()</a>, <a href="#kron">kron()</a>, <a href="#normalise">normalise()</a>, <a href="#repelem">repelem()</a>, <a href="#replace_standalone">replace()</a>, <a href="#repmat">repmat()</a>, <a href="#reshape">reshape()</a>, <a href="#resize">resize()</a>, <a href="#reverse">reverse()</a>, <a href="#shift">shift()</a>, <a href="#symmat">symmatu()/symmatl()</a>, <a href="#trimat">trimatu()/trimatl()</a>, <a href="#t_st_members">.t()</a>, <a href="#trans">trans()</a></li>
+<li>matrix valued functions of matrices: <a href="#clamp">clamp()</a>, <a href="#diagmat">diagmat()</a>, <a href="#diags_spdiags">spdiags()</a>, <a href="#flip">flipud()/fliplr()</a>, <a href="#join">join_rows()</a>, <a href="#join">join_cols()</a>, <a href="#kron">kron()</a>, <a href="#normalise">normalise()</a>, <a href="#repelem">repelem()</a>, <a href="#replace_standalone">replace()</a>, <a href="#repmat">repmat()</a>, <a href="#reshape">reshape()</a>, <a href="#resize">resize()</a>, <a href="#reverse">reverse()</a>, <a href="#circshift">circshift()</a>, <a href="#symmat">symmatu()/symmatl()</a>, <a href="#trimat">trimatu()/trimatl()</a>, <a href="#t_st_members">.t()</a>, <a href="#trans">trans()</a></li>
 <li>generated matrices: <a href="#speye">speye()</a>, <a href="#spones">spones()</a>, <a href="#sprandu_sprandn">sprandu()</a>, <a href="#sprandu_sprandn">sprandn()</a>, <a href="#zeros_standalone">zeros()</a></li>
 <li>eigen decompositions and SVD: <a href="#eigs_sym">eigs_sym()</a>, <a href="#eigs_gen">eigs_gen()</a>, <a href="#svds">svds()</a></li>
 <li>solution of sparse linear systems: <a href="#spsolve">spsolve()</a>
@@ -7994,6 +7967,7 @@ fmat B2 = randu&lt;fmat&gt;(5, 6, distr_param(10,20));
 <ul>
 <li><a href="#randu_randn_member">.randu()</a> (member function)</li>
 <li><a href="#randn">randn()</a></li>
+<li><a href="#rande">rande()</a></li>
 <li><a href="#randg">randg()</a></li>
 <li><a href="#randi">randi()</a></li>
 <li><a href="#imbue">.imbue()</a></li>
@@ -8078,6 +8052,7 @@ fmat B2 = randn&lt;fmat&gt;(5, 6, distr_param(10,5));
 <ul>
 <li><a href="#randu_randn_member">.randn()</a> (member function)</li>
 <li><a href="#randu">randu()</a></li>
+<li><a href="#rande">rande()</a></li>
 <li><a href="#randg">randg()</a></li>
 <li><a href="#randi">randi()</a></li>
 <li><a href="#mvnrnd">mvnrnd()</a></li>
@@ -8092,6 +8067,96 @@ fmat B2 = randn&lt;fmat&gt;(5, 6, distr_param(10,5));
 <br>
 </ul>
 
+<div class="pagebreak"></div><div class="noprint"><hr class="greyline"><br></div>
+<a name="rande"></a>
+<b>rande( )</b> 
+<br><b>rande( distr_param(&thinsp;&lambda;&thinsp;) )</b>
+<br>
+<br><b>rande( n_elem )</b>
+<br><b>rande( n_elem, distr_param(&thinsp;&lambda;&thinsp;) )</b>
+<br>
+<br><b>rande( n_rows, n_cols )</b>
+<br><b>rande( n_rows, n_cols, distr_param(&thinsp;&lambda;&thinsp;) )</b>
+<br>
+<br><b>rande( n_rows, n_cols, n_slices )</b>
+<br><b>rande( n_rows, n_cols, n_slices, distr_param(&thinsp;&lambda;&thinsp;) )</b>
+<br>
+<br><b>rande( size(X) )</b>
+<br><b>rande( size(X), distr_param(&thinsp;&lambda;&thinsp;) )</b>
+<ul>
+<li>
+Generate a scalar, vector, matrix or cube with the elements set to random values from an exponential distribution:
+<br>
+<br>
+<ul>
+<table style="text-align: left;" border="0" cellpadding="0" cellspacing="0">
+<tbody>
+<tr>
+<td>p(&thinsp;x&thinsp;|&thinsp;&lambda;&thinsp;)</td><td>&nbsp;=&nbsp;</td><td>&nbsp;&lambda;&nbsp;exp(&thinsp;-&thinsp;&lambda;&thinsp;x&thinsp;)</td><td>&nbsp;&nbsp;&nbsp;</td><td>for x</td><td>&nbsp;&geq;&nbsp;</td><td>0</td>
+</tr>
+<tr>
+<td>p(&thinsp;x&thinsp;|&thinsp;&lambda;&thinsp;)</td><td>&nbsp;=&nbsp;</td><td>&nbsp;0                                                            </td><td>&nbsp;&nbsp;&nbsp;</td><td>for x</td><td>&nbsp;&lt;&nbsp;</td><td>0</td>
+</tr>
+</tbody>
+</table>
+</ul>
+<br>
+where <i>&lambda;</i> is the rate of change (aka <i>rate parameter</i>), with constraint <i>&lambda;&nbsp;&gt;&nbsp;0</i>
+</li>
+<br>
+<li>The default parameter is <i>&lambda;&nbsp;=&nbsp;1</i>
+</li>
+<br>
+<li>
+Usage:
+<ul>
+<li><i>scalar_type</i> s = rande&lt;<i>scalar_type</i>&gt;( ), where <i>scalar_type</i> is either <i>float</i> or <i>double</i></li>
+<li><i>scalar_type</i> s = rande&lt;<i>scalar_type</i>&gt;( distr_param(&thinsp;&lambda;&thinsp;) ), where <i>scalar_type</i> is either <i>float</i> or <i>double</i></li>
+<br>
+<li><i>vector_type</i> v = rande&lt;<i>vector_type</i>&gt;( n_elem )</li>
+<li><i>vector_type</i> v = rande&lt;<i>vector_type</i>&gt;( n_elem, distr_param(&thinsp;&lambda;&thinsp;) )</li>
+<br>
+<li><i>matrix_type</i> X = rande&lt;<i>matrix_type</i>&gt;( n_rows, n_cols )</li>
+<li><i>matrix_type</i> X = rande&lt;<i>matrix_type</i>&gt;( n_rows, n_cols, distr_param(&thinsp;&lambda;&thinsp;) )</li>
+<br>
+<li><i>cube_type</i> Q = rande&lt;<i>cube_type</i>&gt;( n_rows, n_cols, n_slices )</li>
+<li><i>cube_type</i> Q = rande&lt;<i>cube_type</i>&gt;( n_rows, n_cols, n_slices, distr_param(&thinsp;&lambda;&thinsp;) )</li>
+</ul>
+</li>
+<br>
+<li>
+Examples:
+<ul>
+<pre>
+vec v1 = rande(100);
+vec v2 = rande(100, distr_param(2));
+
+rowvec r1 = rande&lt;rowvec&gt;(100);
+rowvec r2 = rande&lt;rowvec&gt;(100, distr_param(2));
+
+mat A1 = rande(10, 10);
+mat A2 = rande(10, 10, distr_param(2));
+
+fmat B1 = rande&lt;fmat&gt;(10, 10);
+fmat B2 = rande&lt;fmat&gt;(10, 10, distr_param(2));
+</pre>
+</ul>
+</li>
+<li>See also:
+<ul>
+<li><a href="#randg">randg()</a></li>
+<li><a href="#randu">randu()</a></li>
+<li><a href="#randn">randn()</a></li>
+<li><a href="#randi">randi()</a></li>
+<li><a href="#size">size()</a></li>
+<li><a href="#rng_seed">RNG seed setting</a></li>
+<li><a href="https://mathworld.wolfram.com/ExponentialDistribution.html">exponential distribution in MathWorld</a></li>
+<li><a href="https://en.wikipedia.org/wiki/Exponential_distribution">exponential distribution in Wikipedia</a></li>
+</ul>
+</li>
+<br>
+</ul>
+
 <div class="pagebreak"></div><div class="noprint"><hr class="greyline"><br></div>
 <a name="randg"></a>
 <b>randg( )</b> 
@@ -8118,7 +8183,7 @@ Generate a scalar, vector, matrix or cube with the elements set to random values
       <td style="vertical-align: top;">&nbsp;</td><td style="vertical-align: top;">&nbsp;</td><td style="vertical-align: top;"><i>x<sup>&thinsp;a-1</sup> exp(&thinsp;-x&thinsp;/&thinsp;b&thinsp;)</i></td>
     </tr>
     <tr>
-      <td style="vertical-align: top;"><i>p<font size=+1>(</font>x&thinsp;|&thinsp;a,b<font size=+1>)</font></i></td><td style="vertical-align: top;">&nbsp;<font size=+1>=</font>&nbsp;</td><td style="vertical-align: top;"><font size=+1><b><hr></b></font></td>
+      <td style="vertical-align: top;"><i>p<font size=+1>(</font>&thinsp;x&thinsp;|&thinsp;a,b&thinsp;<font size=+1>)</font></i></td><td style="vertical-align: top;">&nbsp;<font size=+1>=</font>&nbsp;</td><td style="vertical-align: top;"><font size=+1><b><hr></b></font></td>
     </tr>
     <tr>
       <td style="vertical-align: top;">&nbsp;</td><td style="vertical-align: top;">&nbsp;</td><td style="vertical-align: top; text-align: center;"><i>b<sup>&thinsp;a</sup> &Gamma;(a)</i></td>
@@ -8169,6 +8234,7 @@ fmat B2 = randg&lt;fmat&gt;(10, 10, distr_param(2,1));
 </li>
 <li>See also:
 <ul>
+<li><a href="#rande">rande()</a></li>
 <li><a href="#randu">randu()</a></li>
 <li><a href="#randn">randn()</a></li>
 <li><a href="#randi">randi()</a></li>
@@ -8176,6 +8242,7 @@ fmat B2 = randg&lt;fmat&gt;(10, 10, distr_param(2,1));
 <li><a href="#imbue">.imbue()</a></li>
 <li><a href="#size">size()</a></li>
 <li><a href="#rng_seed">RNG seed setting</a></li>
+<li><a href="https://mathworld.wolfram.com/GammaDistribution.html">gamma distribution in MathWorld</a></li>
 <li><a href="https://en.wikipedia.org/wiki/Gamma_distribution">gamma distribution in Wikipedia</a></li>
 </ul>
 </li>
@@ -8882,6 +8949,51 @@ See also:
 <br>
 </ul>
 
+<div class="pagebreak"></div><div class="noprint"><hr class="greyline"><br></div>
+<a name="shift"></a>
+<a name="circshift"></a>
+<b>circshift( V, N )</b>
+<br><b>circshift( X, N )</b>
+<br><b>circshift( X, N, dim )</b>
+<ul>
+<li>
+For vector <i>V</i>, generate a copy of the vector with the elements shifted by <i>N</i> positions in a circular manner
+</li>
+<br>
+<li>
+For matrix <i>X</i>, generate a copy of the matrix with the elements shifted by <i>N</i> positions in each column (<i>dim&thinsp;=&thinsp;0</i>), or each row (<i>dim&thinsp;=&thinsp;1</i>)
+</li>
+<br>
+<li>
+<i>N</i> can be positive or negative
+</li>
+<br>
+<li>
+The <i>dim</i> argument is optional; by default <i>dim&thinsp;=&thinsp;0</i> is used
+</li>
+<br>
+<li>
+Examples:
+<ul>
+<pre>
+mat A(4, 5, fill::randu);
+mat B = circshift(A, -1);
+mat C = circshift(A, +1);
+</pre>
+</ul>
+</li>
+<br>
+<li>
+See also:
+<ul>
+<li><a href="#shuffle">shuffle()</a></li>
+<li><a href="#flip">fliplr() &amp; flipud()</a></li>
+<li><a href="#reverse">reverse()</a></li>
+</ul>
+</li>
+<br>
+</ul>
+
 <div class="pagebreak"></div><div class="noprint"><hr class="greyline"><br></div>
 <a name="clamp"></a>
 <b>clamp( X, min_val, max_val )</b>
@@ -9921,7 +10033,7 @@ mat C = flipud(A);
 See also:
 <ul>
 <li><a href="#reverse">reverse()</a></li>
-<li><a href="#shift">shift()</a></li>
+<li><a href="#circshift">circshift()</a></li>
 <li><a href="#swap_rows">.swap_rows() &amp; .swap_cols()</a></li>
 <li><a href="#t_st_members">.t()</a></li>
 </ul>
@@ -11605,7 +11717,7 @@ mat C = reverse(A,1);
 See also:
 <ul>
 <li><a href="#flip">fliplr() &amp; flipud()</a></li>
-<li><a href="#shift">shift()</a></li>
+<li><a href="#circshift">circshift()</a></li>
 <li><a href="#sort">sort()</a></li>
 <li><a href="#trans">trans()</a></li>
 <li><a href="#t_st_members">.t()</a></li>
@@ -11667,50 +11779,6 @@ See also:
 <br>
 </ul>
 
-<div class="pagebreak"></div><div class="noprint"><hr class="greyline"><br></div>
-<a name="shift"></a>
-<b>shift( V, N )</b>
-<br><b>shift( X, N )</b>
-<br><b>shift( X, N, dim )</b>
-<ul>
-<li>
-For vector <i>V</i>, generate a copy of the vector with the elements shifted by <i>N</i> positions in a circular manner
-</li>
-<br>
-<li>
-For matrix <i>X</i>, generate a copy of the matrix with the elements shifted by <i>N</i> positions in each column (<i>dim&thinsp;=&thinsp;0</i>), or each row (<i>dim&thinsp;=&thinsp;1</i>)
-</li>
-<br>
-<li>
-<i>N</i> can be positive or negative
-</li>
-<br>
-<li>
-The <i>dim</i> argument is optional; by default <i>dim&thinsp;=&thinsp;0</i> is used
-</li>
-<br>
-<li>
-Examples:
-<ul>
-<pre>
-mat A(4, 5, fill::randu);
-mat B = shift(A, -1);
-mat C = shift(A, +1);
-</pre>
-</ul>
-</li>
-<br>
-<li>
-See also:
-<ul>
-<li><a href="#shuffle">shuffle()</a></li>
-<li><a href="#flip">fliplr() &amp; flipud()</a></li>
-<li><a href="#reverse">reverse()</a></li>
-</ul>
-</li>
-<br>
-</ul>
-
 <div class="pagebreak"></div><div class="noprint"><hr class="greyline"><br></div>
 <a name="shuffle"></a>
 <b>shuffle( V )</b>
@@ -11742,7 +11810,7 @@ mat B = shuffle(A);
 <li>
 See also:
 <ul>
-<li><a href="#shift">shift()</a></li>
+<li><a href="#circshift">circshift()</a></li>
 <li><a href="#sort">sort()</a></li>
 <li><a href="#unique">unique()</a></li>
 <li><a href="#randi">randi()</a></li>
@@ -11876,9 +11944,6 @@ See also:
 <a name="sort_index"></a>
 <b>sort_index( X )</b>
 <br><b>sort_index( X, sort_direction )</b>
-<br>
-<br><b>stable_sort_index( X )</b>
-<br><b>stable_sort_index( X, sort_direction )</b>
 <ul>
 <li>Return a vector which describes the sorted order of the elements of <i>X</i>
 (ie. it contains the indices of the elements of <i>X</i>)
@@ -11894,8 +11959,6 @@ See also:
 <br>
 <li>The <i>sort_direction</i> argument is optional; <i>sort_direction</i> is either <code>"ascend"</code> or <code>"descend"</code>; by default <code>"ascend"</code> is used</li>
 <br>
-<li>The <i>stable_sort_index()</i> variant preserves the relative order of elements with equivalent values</li>
-<br>
 <li>For matrices and vectors with complex numbers, sorting is via absolute values</li>
 <br>
 <li>
@@ -20748,13 +20811,22 @@ List of additions and changes for each version:
 
 
 
+<a name="version_152"></a>
+<li>Version 15.2:
+<ul>
+<li>added <a href="#rande">rande()</a> for generating matrices with elements from exponential distributions</li>
+<li>shift() has been deprecated in favour of <a href="#circshift">circshift()</a>, for consistency with Matlab/Octave</li>
+<li>reworked detection of aliasing, leading to more efficient compiled code</li>
+</ul>
+</li>
+<br>
 <a name="version_150"></a>
 <li>Version 15.0:
 <ul>
 <li>C++14 is now the minimum required C++ standard</li>
 <li>added preliminary support for matrices with <a href="#fp16_type">half-precision fp16</a> element type</li>
 <li>added second form of <a href="#cond">cond()</a> to allow detection of failures</li>
-<li>added <a href="#repcube">repcube()</a> </li>
+<li>added <a href="#repcube">repcube()</a></li>
 <li>added <code>.freeze()</code> and <code>.unfreeze()</code> member functions to <a href="#wall_clock">wall_clock</a></li>
 <li>extended <a href="#conv">conv()</a> and <a href="#conv2">conv2()</a> to accept the <code>"valid"</code> shape argument</li>
 </ul>
@@ -21448,7 +21520,7 @@ see the <a href="https://arma.sourceforge.net/faq.html#license">Questions page</
 <li>added handling of <a href="#diag">diagonal views</a> by sparse matrices
 <li>expanded <a href="#repmat">repmat()</a> to handle sparse matrices
 <li>expanded <a href="#join">join_rows()</a> and <a href="#join">join_cols()</a> to handle sparse matrices
-<li><a href="#sort_index">sort_index()</a> and <a href="#sort_index">stable_sort_index()</a> have been placed in the delayed operations framework for increased efficiency
+<li><a href="#sort_index">sort_index()</a> has been placed in the delayed operations framework for increased efficiency
 <li>use of <a href="#uword">64 bit integers</a> is automatically enabled when using a C++11 compiler</li>
 </ul>
 </li>
@@ -21660,7 +21732,7 @@ eg. <i>mat&nbsp;X(4,&nbsp;5,&nbsp;fill::zeros)</i></li>
 <li>faster handling of compound expressions with submatrices and subcubes</li>
 <li>faster <a href="#trace">trace()</a></li>
 <li>added support for loading matrices as text files with <i>NaN</i> and <i>Inf</i> elements</li>
-<li>added <a href="#sort_index">stable_sort_index()</a>, which preserves the relative order of elements with equivalent values</li>
+<!-- <li>added <a href="#sort_index">stable_sort_index()</a>, which preserves the relative order of elements with equivalent values</li> -->
 <li>added handling of <a href="#SpMat">sparse matrices</a> by <a href="#stats_fns">mean()</a>, <a href="#stats_fns">var()</a>, <a href="#norm">norm()</a>, <a href="#abs">abs()</a>, <a href="#misc_fns">square()</a>, <a href="#misc_fns">sqrt()</a></li>
 <li>added saving and loading of sparse matrices in <i>arma_binary</i> format</li>
 </ul>

include/armadillo

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

include/armadillo_bits/Base_meat.hpp

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

include/armadillo_bits/Col_bones.hpp

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

include/armadillo_bits/Col_meat.hpp

@@ -483,7 +483,7 @@ Col<eT>::operator=(Col<eT>&& X)
 //   if(X.n_cols != 1)  { const Mat<eT>& XX = X; Mat<eT>::operator=(XX); return *this; }
 //   
 //   (*this).steal_mem(X, true);
-//   
+//    
 //   return *this;
 //   }
 
@@ -521,11 +521,9 @@ template<typename eT>
 template<typename T1>
 inline
 Col<eT>::Col(const Base<eT,T1>& X)
-  : Mat<eT>(arma_vec_indicator(), 1)
+  : Mat<eT>(X.get_ref(), arma_vec_indicator(), 1)
   {
-  arma_debug_sigprint();
-  
-  Mat<eT>::operator=(X.get_ref());
+  arma_debug_sigprint_this(this);
   }
 
 
@@ -549,11 +547,9 @@ template<typename eT>
 template<typename T1>
 inline
 Col<eT>::Col(const SpBase<eT,T1>& X)
-  : Mat<eT>(arma_vec_indicator(), 1)
+  : Mat<eT>(X.get_ref(), arma_vec_indicator(), 1)
   {
   arma_debug_sigprint_this(this);
-  
-  Mat<eT>::operator=(X.get_ref());
   }
 
 
@@ -607,12 +603,9 @@ Col<eT>::Col
   const Base<typename Col<eT>::pod_type, T1>& A,
   const Base<typename Col<eT>::pod_type, T2>& B
   )
+  : Mat<eT>(A.get_ref(), B.get_ref(), arma_vec_indicator(), 1)
   {
-  arma_debug_sigprint();
-  
-  access::rw(Mat<eT>::vec_state) = 1;
-  
-  Mat<eT>::init(A,B);
+  arma_debug_sigprint_this(this);
   }
 
 
@@ -621,12 +614,9 @@ template<typename eT>
 template<typename T1>
 inline
 Col<eT>::Col(const BaseCube<eT,T1>& X)
+  : Mat<eT>(X.get_ref(), arma_vec_indicator(), 1)
   {
-  arma_debug_sigprint();
-  
-  access::rw(Mat<eT>::vec_state) = 1;
-  
-  Mat<eT>::operator=(X);
+  arma_debug_sigprint_this(this);
   }
 
 
@@ -649,12 +639,9 @@ Col<eT>::operator=(const BaseCube<eT,T1>& X)
 template<typename eT>
 inline
 Col<eT>::Col(const subview_cube<eT>& X)
+  : Mat<eT>(X, arma_vec_indicator(), 1)
   {
-  arma_debug_sigprint();
-  
-  access::rw(Mat<eT>::vec_state) = 1;
-  
-  Mat<eT>::operator=(X);
+  arma_debug_sigprint_this(this);
   }
 
 
@@ -1277,6 +1264,13 @@ Col<eT>::end_row(const uword row_num) const
 
 
 
+
+
+
+//
+
+
+
 template<typename eT>
 template<uword fixed_n_elem>
 arma_inline
@@ -1872,6 +1866,10 @@ Col<eT>::fixed<fixed_n_elem>::ones()
 
 
 
+//
+
+
+
 template<typename eT>
 inline
 Col<eT>::Col(const arma_fixed_indicator&, const uword in_n_elem, const eT* in_mem)

include/armadillo_bits/Cube_meat.hpp

@@ -3148,7 +3148,7 @@ Cube<eT>::Cube(const GlueCube<T1, T2, glue_type>& X)
   {
   arma_debug_sigprint_this(this);
   
-  this->operator=(X);
+  (*this).operator=(X);
   }
 
 

include/armadillo_bits/Mat_bones.hpp

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

include/armadillo_bits/Mat_meat.hpp

@@ -880,7 +880,26 @@ Mat<eT>::operator/=(const eT val)
 
 
 
-//! construct a matrix from a given matrix
+template<typename eT>
+inline
+Mat<eT>::Mat(const Mat<eT>& in_mat, const arma_vec_indicator&, const uhword in_vec_state)
+  : n_rows( (in_vec_state == 2) ? 1 : 0 )
+  , n_cols( (in_vec_state == 1) ? 1 : 0 )
+  , n_elem(0)
+  , n_alloc(0)
+  , vec_state(in_vec_state)
+  , mem_state(0)
+  , mem()
+  {
+  arma_debug_sigprint(arma_str::format("this: %x; in_mat: %x") % this % &in_mat);
+  
+  init_warm(in_mat.n_rows, in_mat.n_cols);
+  
+  arrayops::copy( memptr(), in_mat.mem, in_mat.n_elem );
+  }
+
+
+
 template<typename eT>
 inline
 Mat<eT>::Mat(const Mat<eT>& in_mat)
@@ -901,7 +920,6 @@ Mat<eT>::Mat(const Mat<eT>& in_mat)
 
 
 
-//! construct a matrix from a given matrix
 template<typename eT>
 inline
 Mat<eT>&
@@ -1496,6 +1514,25 @@ Mat<eT>::operator/=(const Mat<eT>& m)
 
 
 
+template<typename eT>
+template<typename T1>
+inline
+Mat<eT>::Mat(const BaseCube<eT,T1>& X, const arma_vec_indicator&, const uhword in_vec_state)
+  : n_rows( (in_vec_state == 2) ? 1 : 0 )
+  , n_cols( (in_vec_state == 1) ? 1 : 0 )
+  , n_elem(0)
+  , n_alloc(0)
+  , vec_state(in_vec_state)
+  , mem_state(0)
+  , mem()
+  {
+  arma_debug_sigprint_this(this);
+  
+  (*this).operator=(X);
+  }
+
+
+
 template<typename eT>
 template<typename T1>
 inline
@@ -1943,15 +1980,30 @@ Mat<eT>::operator/=(const BaseCube<eT,T1>& X)
 
 
 
+template<typename eT>
+template<typename T1, typename T2>
+inline
+Mat<eT>::Mat(const Base<typename Mat<eT>::pod_type,T1>& A, const Base<typename Mat<eT>::pod_type,T2>& B, const arma_vec_indicator&, const uhword in_vec_state)
+  : n_rows( (in_vec_state == 2) ? 1 : 0 )
+  , n_cols( (in_vec_state == 1) ? 1 : 0 )
+  , n_elem(0)
+  , n_alloc(0)
+  , vec_state(in_vec_state)
+  , mem_state(0)
+  , mem()
+  {
+  arma_debug_sigprint_this(this);
+  
+  init(A,B);
+  }
+
+
+
 //! for constructing a complex matrix out of two non-complex matrices
 template<typename eT>
 template<typename T1, typename T2>
 inline
-Mat<eT>::Mat
-  (
-  const Base<typename Mat<eT>::pod_type,T1>& A,
-  const Base<typename Mat<eT>::pod_type,T2>& B
-  )
+Mat<eT>::Mat(const Base<typename Mat<eT>::pod_type,T1>& A, const Base<typename Mat<eT>::pod_type,T2>& B)
   : n_rows(0)
   , n_cols(0)
   , n_elem(0)
@@ -1969,18 +2021,18 @@ Mat<eT>::Mat
 
 template<typename eT>
 inline
-Mat<eT>::Mat(const subview<eT>& X, const bool use_colmem)
+Mat<eT>::Mat(const subview<eT>& X, const bool reuse_mem)
   : n_rows(X.n_rows)
   , n_cols(X.n_cols)
   , n_elem(X.n_elem)
   , n_alloc(0)
   , vec_state(0)
-  , mem_state(use_colmem ? 3 : 0)
-  , mem      (use_colmem ? X.colptr(0) : nullptr)
+  , mem_state(reuse_mem ? 3           : 0      )
+  , mem      (reuse_mem ? X.colptr(0) : nullptr)
   {
   arma_debug_sigprint_this(this);
   
-  if(use_colmem)
+  if(reuse_mem)
     {
     arma_debug_print("Mat::Mat(): using existing memory in a submatrix");
     }
@@ -1994,6 +2046,26 @@ Mat<eT>::Mat(const subview<eT>& X, const bool use_colmem)
 
 
 
+template<typename eT>
+inline
+Mat<eT>::Mat(const subview<eT>& X, const arma_vec_indicator&, const uhword in_vec_state)
+  : n_rows( (in_vec_state == 2) ? 1 : 0 )
+  , n_cols( (in_vec_state == 1) ? 1 : 0 )
+  , n_elem(0)
+  , n_alloc(0)
+  , vec_state(in_vec_state)
+  , mem_state(0)
+  , mem()
+  {
+  arma_debug_sigprint_this(this);
+  
+  init_warm(X.n_rows, X.n_cols);
+    
+  subview<eT>::extract(*this, X);
+  }
+
+
+
 //! construct a matrix from subview (eg. construct a matrix from a delayed submatrix operation)
 template<typename eT>
 inline
@@ -2197,6 +2269,24 @@ Mat<eT>::Mat(const xtrans_mat<eT,do_conj>& X)
 
 
 
+template<typename eT>
+inline
+Mat<eT>::Mat(const subview_cube<eT>& x, const arma_vec_indicator&, const uhword in_vec_state)
+  : n_rows( (in_vec_state == 2) ? 1 : 0 )
+  , n_cols( (in_vec_state == 1) ? 1 : 0 )
+  , n_elem(0)
+  , n_alloc(0)
+  , vec_state(in_vec_state)
+  , mem_state(0)
+  , mem()
+  {
+  arma_debug_sigprint_this(this);
+  
+  (*this).operator=(x);
+  }
+
+
+
 //! construct a matrix from a subview_cube instance
 template<typename eT>
 inline
@@ -2211,7 +2301,7 @@ Mat<eT>::Mat(const subview_cube<eT>& x)
   {
   arma_debug_sigprint_this(this);
   
-  this->operator=(x);
+  (*this).operator=(x);
   }
 
 
@@ -2308,6 +2398,26 @@ Mat<eT>::operator/=(const subview_cube<eT>& X)
 
 
 
+template<typename eT>
+inline
+Mat<eT>::Mat(const diagview<eT>& X, const arma_vec_indicator&, const uhword in_vec_state)
+  : n_rows( (in_vec_state == 2) ? 1 : 0 )
+  , n_cols( (in_vec_state == 1) ? 1 : 0 )
+  , n_elem(0)
+  , n_alloc(0)
+  , vec_state(in_vec_state)
+  , mem_state(0)
+  , mem()
+  {
+  arma_debug_sigprint_this(this);
+  
+  init_warm(X.n_rows, X.n_cols);
+  
+  diagview<eT>::extract(*this, X);
+  }
+
+
+
 //! construct a matrix from diagview (eg. construct a matrix from a delayed diag operation)
 template<typename eT>
 inline
@@ -2432,6 +2542,25 @@ Mat<eT>::operator/=(const diagview<eT>& X)
 
 
 
+template<typename eT>
+template<typename T1>
+inline
+Mat<eT>::Mat(const subview_elem1<eT,T1>& X, const arma_vec_indicator&, const uhword in_vec_state)
+  : n_rows( (in_vec_state == 2) ? 1 : 0 )
+  , n_cols( (in_vec_state == 1) ? 1 : 0 )
+  , n_elem(0)
+  , n_alloc(0)
+  , vec_state(in_vec_state)
+  , mem_state(0)
+  , mem()
+  {
+  arma_debug_sigprint_this(this);
+  
+  (*this).operator=(X);
+  }
+
+
+
 template<typename eT>
 template<typename T1>
 inline
@@ -2446,7 +2575,7 @@ Mat<eT>::Mat(const subview_elem1<eT,T1>& X)
   {
   arma_debug_sigprint_this(this);
   
-  this->operator=(X);
+  (*this).operator=(X);
   }
 
 
@@ -2541,6 +2670,25 @@ Mat<eT>::operator/=(const subview_elem1<eT,T1>& X)
 
 
 
+template<typename eT>
+template<typename T1, typename T2>
+inline
+Mat<eT>::Mat(const subview_elem2<eT,T1,T2>& X, const arma_vec_indicator&, const uhword in_vec_state)
+  : n_rows( (in_vec_state == 2) ? 1 : 0 )
+  , n_cols( (in_vec_state == 1) ? 1 : 0 )
+  , n_elem(0)
+  , n_alloc(0)
+  , vec_state(in_vec_state)
+  , mem_state(0)
+  , mem()
+  {
+  arma_debug_sigprint_this(this);
+  
+  (*this).operator=(X);
+  }
+
+
+
 template<typename eT>
 template<typename T1, typename T2>
 inline
@@ -2555,7 +2703,7 @@ Mat<eT>::Mat(const subview_elem2<eT,T1,T2>& X)
   {
   arma_debug_sigprint_this(this);
   
-  this->operator=(X);
+  (*this).operator=(X);
   }
 
 
@@ -2650,6 +2798,25 @@ Mat<eT>::operator/=(const subview_elem2<eT,T1,T2>& X)
 
 
 
+template<typename eT>
+template<typename T1>
+inline
+Mat<eT>::Mat(const SpBase<eT, T1>& m, const arma_vec_indicator&, const uhword in_vec_state)
+  : n_rows( (in_vec_state == 2) ? 1 : 0 )
+  , n_cols( (in_vec_state == 1) ? 1 : 0 )
+  , n_elem(0)
+  , n_alloc(0)
+  , vec_state(in_vec_state)
+  , mem_state(0)
+  , mem()
+  {
+  arma_debug_sigprint_this(this);
+  
+  (*this).operator=(m);
+  }
+
+
+
 template<typename eT>
 template<typename T1>
 inline
@@ -2835,6 +3002,24 @@ Mat<eT>::operator/=(const SpBase<eT, T1>& m)
 
 
 
+template<typename eT>
+inline
+Mat<eT>::Mat(const SpSubview<eT>& X, const arma_vec_indicator&, const uhword in_vec_state)
+  : n_rows( (in_vec_state == 2) ? 1 : 0 )
+  , n_cols( (in_vec_state == 1) ? 1 : 0 )
+  , n_elem(0)
+  , n_alloc(0)
+  , vec_state(in_vec_state)
+  , mem_state(0)
+  , mem()
+  {
+  arma_debug_sigprint_this(this);
+  
+  (*this).operator=(X);
+  }
+
+
+
 template<typename eT>
 inline
 Mat<eT>::Mat(const SpSubview<eT>& X)
@@ -3018,6 +3203,26 @@ Mat<eT>::operator-=(const SpSubview<eT>& X)
 
 
 
+template<typename eT>
+inline
+Mat<eT>::Mat(const spdiagview<eT>& X, const arma_vec_indicator&, const uhword in_vec_state)
+  : n_rows( (in_vec_state == 2) ? 1 : 0 )
+  , n_cols( (in_vec_state == 1) ? 1 : 0 )
+  , n_elem(0)
+  , n_alloc(0)
+  , vec_state(in_vec_state)
+  , mem_state(0)
+  , mem()
+  {
+  arma_debug_sigprint_this(this);
+  
+  init_warm(X.n_rows, X.n_cols);
+  
+  spdiagview<eT>::extract(*this, X);
+  }
+
+
+
 template<typename eT>
 inline
 Mat<eT>::Mat(const spdiagview<eT>& X)
@@ -4925,6 +5130,29 @@ Mat<eT>::insert_cols(const uword col_num, const Base<eT,T1>& X)
 
 
 
+template<typename eT>
+template<typename T1, typename gen_type>
+inline
+Mat<eT>::Mat(const Gen<T1, gen_type>& X, const arma_vec_indicator&, const uhword in_vec_state)
+  : n_rows( (in_vec_state == 2) ? 1 : 0 )
+  , n_cols( (in_vec_state == 1) ? 1 : 0 )
+  , n_elem(0)
+  , n_alloc(0)
+  , vec_state(in_vec_state)
+  , mem_state(0)
+  , mem()
+  {
+  arma_debug_sigprint_this(this);
+  
+  arma_type_check(( is_same_type< eT, typename T1::elem_type >::no ));
+  
+  init_warm(X.n_rows, X.n_cols);
+  
+  X.apply(*this);
+  }
+
+
+
 template<typename eT>
 template<typename T1, typename gen_type>
 inline
@@ -5052,6 +5280,27 @@ Mat<eT>::operator/=(const Gen<T1, gen_type>& X)
 
 
 
+template<typename eT>
+template<typename T1, typename op_type>
+inline
+Mat<eT>::Mat(const Op<T1, op_type>& X, const arma_vec_indicator&, const uhword in_vec_state)
+  : n_rows( (in_vec_state == 2) ? 1 : 0 )
+  , n_cols( (in_vec_state == 1) ? 1 : 0 )
+  , n_elem(0)
+  , n_alloc(0)
+  , vec_state(in_vec_state)
+  , mem_state(0)
+  , mem()
+  {
+  arma_debug_sigprint_this(this);
+  
+  arma_type_check(( is_same_type< eT, typename T1::elem_type >::no ));
+  
+  op_type::apply(static_cast< Mat_noalias<eT>& >(*this), X);
+  }
+
+
+
 //! create a matrix from Op, ie. run the previously delayed unary operations
 template<typename eT>
 template<typename T1, typename op_type>
@@ -5066,10 +5315,10 @@ Mat<eT>::Mat(const Op<T1, op_type>& X)
   , mem()
   {
   arma_debug_sigprint_this(this);
-
+  
   arma_type_check(( is_same_type< eT, typename T1::elem_type >::no ));
   
-  op_type::apply(*this, X);
+  op_type::apply(static_cast< Mat_noalias<eT>& >(*this), X);
   }
 
 
@@ -5082,7 +5331,7 @@ Mat<eT>&
 Mat<eT>::operator=(const Op<T1, op_type>& X)
   {
   arma_debug_sigprint();
-
+  
   arma_type_check(( is_same_type< eT, typename T1::elem_type >::no ));
   
   op_type::apply(*this, X);
@@ -5182,6 +5431,37 @@ Mat<eT>::operator/=(const Op<T1, op_type>& X)
 
 
 
+template<typename eT>
+template<typename T1, typename eop_type>
+inline
+Mat<eT>::Mat(const eOp<T1, eop_type>& X, const arma_vec_indicator&, const uhword in_vec_state)
+  : n_rows( (in_vec_state == 2) ? 1 : 0 )
+  , n_cols( (in_vec_state == 1) ? 1 : 0 )
+  , n_elem(0)
+  , n_alloc(0)
+  , vec_state(in_vec_state)
+  , mem_state(0)
+  , mem()
+  {
+  arma_debug_sigprint_this(this);
+  
+  arma_type_check(( is_same_type< eT, typename T1::elem_type >::no ));
+  
+  init_warm(X.get_n_rows(), X.get_n_cols());
+  
+  if(arma_config::optimise_powexpr && is_same_type<eop_type, eop_pow>::value)
+    {
+    constexpr bool eT_ok = is_real_or_cx<eT>::value;
+    
+    if(          X.aux == eT(2)   )  { eop_square::apply(*this, reinterpret_cast< const eOp<T1, eop_square>& >(X)); return; }
+    if(eT_ok && (X.aux == eT(0.5)))  {   eop_sqrt::apply(*this, reinterpret_cast< const eOp<T1, eop_sqrt  >& >(X)); return; }
+    }
+  
+  eop_type::apply(*this, X);
+  }
+
+
+
 //! create a matrix from eOp, ie. run the previously delayed unary operations
 template<typename eT>
 template<typename T1, typename eop_type>
@@ -5379,6 +5659,25 @@ Mat<eT>::operator/=(const eOp<T1, eop_type>& X)
 
 
 
+template<typename eT>
+template<typename T1, typename op_type>
+inline
+Mat<eT>::Mat(const mtOp<eT, T1, op_type>& X, const arma_vec_indicator&, const uhword in_vec_state)
+  : n_rows( (in_vec_state == 2) ? 1 : 0 )
+  , n_cols( (in_vec_state == 1) ? 1 : 0 )
+  , n_elem(0)
+  , n_alloc(0)
+  , vec_state(in_vec_state)
+  , mem_state(0)
+  , mem()
+  {
+  arma_debug_sigprint_this(this);
+  
+  op_type::apply(static_cast< Mat_noalias<eT>& >(*this), X);
+  }
+
+
+
 template<typename eT>
 template<typename T1, typename op_type>
 inline
@@ -5393,7 +5692,7 @@ Mat<eT>::Mat(const mtOp<eT, T1, op_type>& X)
   {
   arma_debug_sigprint_this(this);
   
-  op_type::apply(*this, X);
+  op_type::apply(static_cast< Mat_noalias<eT>& >(*this), X);
   }
 
 
@@ -5488,6 +5787,27 @@ Mat<eT>::operator/=(const mtOp<eT, T1, op_type>& X)
 
 
 
+template<typename eT>
+template<typename T1, typename op_type>
+inline
+Mat<eT>::Mat(const CubeToMatOp<T1, op_type>& X, const arma_vec_indicator&, const uhword in_vec_state)
+  : n_rows( (in_vec_state == 2) ? 1 : 0 )
+  , n_cols( (in_vec_state == 1) ? 1 : 0 )
+  , n_elem(0)
+  , n_alloc(0)
+  , vec_state(in_vec_state)
+  , mem_state(0)
+  , mem()
+  {
+  arma_debug_sigprint_this(this);
+  
+  arma_type_check(( is_same_type< eT, typename T1::elem_type >::no ));
+  
+  op_type::apply(*this, X);
+  }
+
+
+
 template<typename eT>
 template<typename T1, typename op_type>
 inline
@@ -5501,9 +5821,9 @@ Mat<eT>::Mat(const CubeToMatOp<T1, op_type>& X)
   , mem()
   {
   arma_debug_sigprint_this(this);
-
+  
   arma_type_check(( is_same_type< eT, typename T1::elem_type >::no ));
-
+  
   op_type::apply(*this, X);
   }
 
@@ -5516,7 +5836,7 @@ Mat<eT>&
 Mat<eT>::operator=(const CubeToMatOp<T1, op_type>& X)
   {
   arma_debug_sigprint();
-
+  
   arma_type_check(( is_same_type< eT, typename T1::elem_type >::no ));
   
   op_type::apply(*this, X);
@@ -5611,6 +5931,27 @@ Mat<eT>::operator/=(const CubeToMatOp<T1, op_type>& X)
 
 
 
+template<typename eT>
+template<typename T1, typename op_type>
+inline
+Mat<eT>::Mat(const SpToDOp<T1, op_type>& X, const arma_vec_indicator&, const uhword in_vec_state)
+  : n_rows( (in_vec_state == 2) ? 1 : 0 )
+  , n_cols( (in_vec_state == 1) ? 1 : 0 )
+  , n_elem(0)
+  , n_alloc(0)
+  , vec_state(in_vec_state)
+  , mem_state(0)
+  , mem()
+  {
+  arma_debug_sigprint_this(this);
+  
+  arma_type_check(( is_same_type< eT, typename T1::elem_type >::no ));
+  
+  op_type::apply(*this, X);
+  }
+
+
+
 template<typename eT>
 template<typename T1, typename op_type>
 inline
@@ -5624,9 +5965,9 @@ Mat<eT>::Mat(const SpToDOp<T1, op_type>& X)
   , mem()
   {
   arma_debug_sigprint_this(this);
-
+  
   arma_type_check(( is_same_type< eT, typename T1::elem_type >::no ));
-
+  
   op_type::apply(*this, X);
   }
 
@@ -5640,7 +5981,7 @@ Mat<eT>&
 Mat<eT>::operator=(const SpToDOp<T1, op_type>& X)
   {
   arma_debug_sigprint();
-
+  
   arma_type_check(( is_same_type< eT, typename T1::elem_type >::no ));
   
   op_type::apply(*this, X);
@@ -5740,6 +6081,25 @@ Mat<eT>::operator/=(const SpToDOp<T1, op_type>& X)
 
 
 
+template<typename eT>
+template<typename T1, typename op_type>
+inline
+Mat<eT>::Mat(const mtSpReduceOp<eT, T1, op_type>& X, const arma_vec_indicator&, const uhword in_vec_state)
+  : n_rows( (in_vec_state == 2) ? 1 : 0 )
+  , n_cols( (in_vec_state == 1) ? 1 : 0 )
+  , n_elem(0)
+  , n_alloc(0)
+  , vec_state(in_vec_state)
+  , mem_state(0)
+  , mem()
+  {
+  arma_debug_sigprint_this(this);
+  
+  op_type::apply(*this, X);
+  }
+
+
+
 template<typename eT>
 template<typename T1, typename op_type>
 inline
@@ -5753,7 +6113,7 @@ Mat<eT>::Mat(const mtSpReduceOp<eT, T1, op_type>& X)
   , mem()
   {
   arma_debug_sigprint_this(this);
-
+  
   op_type::apply(*this, X);
   }
 
@@ -5766,7 +6126,7 @@ Mat<eT>&
 Mat<eT>::operator=(const mtSpReduceOp<eT, T1, op_type>& X)
   {
   arma_debug_sigprint();
-
+  
   op_type::apply(*this, X);
   
   return *this;
@@ -5849,6 +6209,28 @@ Mat<eT>::operator/=(const mtSpReduceOp<eT, T1, op_type>& X)
 
 
 
+template<typename eT>
+template<typename T1, typename T2, typename glue_type>
+inline
+Mat<eT>::Mat(const Glue<T1, T2, glue_type>& X, const arma_vec_indicator&, const uhword in_vec_state)
+  : n_rows( (in_vec_state == 2) ? 1 : 0 )
+  , n_cols( (in_vec_state == 1) ? 1 : 0 )
+  , n_elem(0)
+  , n_alloc(0)
+  , vec_state(in_vec_state)
+  , mem_state(0)
+  , mem()
+  {
+  arma_debug_sigprint_this(this);
+  
+  arma_type_check(( is_same_type< eT, typename T1::elem_type >::no ));
+  arma_type_check(( is_same_type< eT, typename T2::elem_type >::no ));
+  
+  glue_type::apply(static_cast< Mat_noalias<eT>& >(*this), X);
+  }
+
+
+
 //! create a matrix from Glue, ie. run the previously delayed binary operations
 template<typename eT>
 template<typename T1, typename T2, typename glue_type>
@@ -5867,7 +6249,7 @@ Mat<eT>::Mat(const Glue<T1, T2, glue_type>& X)
   arma_type_check(( is_same_type< eT, typename T1::elem_type >::no ));
   arma_type_check(( is_same_type< eT, typename T2::elem_type >::no ));
   
-  glue_type::apply(*this, X);
+  glue_type::apply(static_cast< Mat_noalias<eT>& >(*this), X);
   }
 
 
@@ -6016,6 +6398,30 @@ Mat<eT>::operator-=(const Glue<T1, T2, glue_times>& X)
 
 
 
+template<typename eT>
+template<typename T1, typename T2, typename eglue_type>
+inline
+Mat<eT>::Mat(const eGlue<T1, T2, eglue_type>& X, const arma_vec_indicator&, const uhword in_vec_state)
+  : n_rows( (in_vec_state == 2) ? 1 : 0 )
+  , n_cols( (in_vec_state == 1) ? 1 : 0 )
+  , n_elem(0)
+  , n_alloc(0)
+  , vec_state(in_vec_state)
+  , mem_state(0)
+  , mem()
+  {
+  arma_debug_sigprint_this(this);
+  
+  arma_type_check(( is_same_type< eT, typename T1::elem_type >::no ));
+  arma_type_check(( is_same_type< eT, typename T2::elem_type >::no ));
+  
+  init_warm(X.get_n_rows(), X.get_n_cols());
+  
+  eglue_type::apply(*this, X);
+  }
+
+
+
 //! create a matrix from eGlue, ie. run the previously delayed binary operations
 template<typename eT>
 template<typename T1, typename T2, typename eglue_type>
@@ -6199,6 +6605,25 @@ Mat<eT>::operator/=(const eGlue<T1, T2, eglue_type>& X)
 
 
 
+template<typename eT>
+template<typename T1, typename T2, typename glue_type>
+inline
+Mat<eT>::Mat(const mtGlue<eT, T1, T2, glue_type>& X, const arma_vec_indicator&, const uhword in_vec_state)
+  : n_rows( (in_vec_state == 2) ? 1 : 0 )
+  , n_cols( (in_vec_state == 1) ? 1 : 0 )
+  , n_elem(0)
+  , n_alloc(0)
+  , vec_state(in_vec_state)
+  , mem_state(0)
+  , mem()
+  {
+  arma_debug_sigprint_this(this);
+  
+  glue_type::apply(static_cast< Mat_noalias<eT>& >(*this), X);
+  }
+
+
+
 template<typename eT>
 template<typename T1, typename T2, typename glue_type>
 inline
@@ -6213,7 +6638,7 @@ Mat<eT>::Mat(const mtGlue<eT, T1, T2, glue_type>& X)
   {
   arma_debug_sigprint_this(this);
   
-  glue_type::apply(*this, X);
+  glue_type::apply(static_cast< Mat_noalias<eT>& >(*this), X);
   }
 
 
@@ -6310,6 +6735,28 @@ Mat<eT>::operator/=(const mtGlue<eT, T1, T2, glue_type>& X)
 
 
 
+template<typename eT>
+template<typename T1, typename T2, typename glue_type>
+inline
+Mat<eT>::Mat(const SpToDGlue<T1, T2, glue_type>& X, const arma_vec_indicator&, const uhword in_vec_state)
+  : n_rows( (in_vec_state == 2) ? 1 : 0 )
+  , n_cols( (in_vec_state == 1) ? 1 : 0 )
+  , n_elem(0)
+  , n_alloc(0)
+  , vec_state(in_vec_state)
+  , mem_state(0)
+  , mem()
+  {
+  arma_debug_sigprint_this(this);
+  
+  arma_type_check(( is_same_type< eT, typename T1::elem_type >::no ));
+  arma_type_check(( is_same_type< eT, typename T2::elem_type >::no ));
+  
+  glue_type::apply(*this, X);
+  }
+
+
+
 template<typename eT>
 template<typename T1, typename T2, typename glue_type>
 inline

include/armadillo_bits/Row_bones.hpp

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

include/armadillo_bits/Row_meat.hpp

@@ -521,11 +521,9 @@ template<typename eT>
 template<typename T1>
 inline
 Row<eT>::Row(const Base<eT,T1>& X)
-  : Mat<eT>(arma_vec_indicator(), 2)
+  : Mat<eT>(X.get_ref(), arma_vec_indicator(), 2)
   {
-  arma_debug_sigprint();
-  
-  Mat<eT>::operator=(X.get_ref());
+  arma_debug_sigprint_this(this);
   }
 
 
@@ -549,11 +547,9 @@ template<typename eT>
 template<typename T1>
 inline
 Row<eT>::Row(const SpBase<eT,T1>& X)
-  : Mat<eT>(arma_vec_indicator(), 2)
+  : Mat<eT>(X.get_ref(), arma_vec_indicator(), 2)
   {
-  arma_debug_sigprint();
-  
-  Mat<eT>::operator=(X.get_ref());
+  arma_debug_sigprint_this(this);
   }
 
 
@@ -607,12 +603,9 @@ Row<eT>::Row
   const Base<typename Row<eT>::pod_type, T1>& A,
   const Base<typename Row<eT>::pod_type, T2>& B
   )
+  : Mat<eT>(A.get_ref(), B.get_ref(), arma_vec_indicator(), 2)
   {
-  arma_debug_sigprint();
-  
-  access::rw(Mat<eT>::vec_state) = 2;
-  
-  Mat<eT>::init(A,B);
+  arma_debug_sigprint_this(this);
   }
 
 
@@ -621,12 +614,9 @@ template<typename eT>
 template<typename T1>
 inline
 Row<eT>::Row(const BaseCube<eT,T1>& X)
+  : Mat<eT>(X.get_ref(), arma_vec_indicator(), 2)
   {
-  arma_debug_sigprint();
-  
-  access::rw(Mat<eT>::vec_state) = 2;
-  
-  Mat<eT>::operator=(X);
+  arma_debug_sigprint_this(this);
   }
 
 
@@ -649,12 +639,9 @@ Row<eT>::operator=(const BaseCube<eT,T1>& X)
 template<typename eT>
 inline
 Row<eT>::Row(const subview_cube<eT>& X)
+  : Mat<eT>(X, arma_vec_indicator(), 2)
   {
-  arma_debug_sigprint();
-  
-  access::rw(Mat<eT>::vec_state) = 2;
-  
-  Mat<eT>::operator=(X);
+  arma_debug_sigprint_this(this);
   }
 
 
@@ -1277,6 +1264,20 @@ Row<eT>::end_row(const uword row_num) const
 
 
 
+template<typename eT>
+inline
+Row<eT>::Row(const subview<eT>& X, const bool reuse_mem)
+  : Mat<eT>(X, reuse_mem)
+  {
+  arma_debug_sigprint_this(this);
+  }
+
+
+
+//
+
+
+
 template<typename eT>
 template<uword fixed_n_elem>
 arma_inline
@@ -1872,6 +1873,10 @@ Row<eT>::fixed<fixed_n_elem>::ones()
 
 
 
+//
+
+
+
 template<typename eT>
 inline
 Row<eT>::Row(const arma_fixed_indicator&, const uword in_n_elem, const eT* in_mem)

include/armadillo_bits/SpMat_bones.hpp

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

include/armadillo_bits/SpMat_meat.hpp

@@ -1666,7 +1666,7 @@ SpMat<eT>::SpMat(const SpOp<T1, spop_type>& X)
   
   arma_type_check(( is_same_type< eT, typename T1::elem_type >::no ));
   
-  spop_type::apply(*this, X);
+  spop_type::apply(static_cast< SpMat_noalias<eT>& >(*this), X);
   
   sync_csc();          // in case apply() used element accessors
   invalidate_cache();  // in case apply() modified the CSC representation
@@ -1806,7 +1806,7 @@ SpMat<eT>::SpMat(const SpGlue<T1, T2, spglue_type>& X)
   
   arma_type_check(( is_same_type< eT, typename T1::elem_type >::no ));
   
-  spglue_type::apply(*this, X);
+  spglue_type::apply(static_cast< SpMat_noalias<eT>& >(*this), X);
   
   sync_csc();          // in case apply() used element accessors
   invalidate_cache();  // in case apply() modified the CSC representation

include/armadillo_bits/arma_forward.hpp

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

include/armadillo_bits/arma_rng.hpp

@@ -116,6 +116,7 @@ struct arma_rng
   template<typename eT> struct randu;
   template<typename eT> struct randn;
   template<typename eT> struct randg;
+  template<typename eT> struct rande;
   };
 
 
@@ -1029,4 +1030,45 @@ struct arma_rng::randg
 
 
 
+//
+
+
+
+template<typename eT>
+struct arma_rng::rande
+  {
+  inline
+  static
+  void
+  fill(eT* mem, const uword N, const double lambda)
+    {
+    #if defined(ARMA_USE_CXX11_RNG)
+      {
+      std::exponential_distribution<double> local_e_distr(lambda);
+      
+      std::mt19937_64& producer = arma_rng::get_producer();
+      
+      arma_rng::lock_producer();
+      
+      for(uword i=0; i<N; ++i)  { mem[i] = eT(local_e_distr(producer)); }
+      
+      arma_rng::unlock_producer();
+      }
+    #else
+      {
+      typedef typename std::mt19937_64::result_type local_seed_type;
+      
+      std::mt19937_64                       local_engine;
+      std::exponential_distribution<double> local_e_distr(lambda);
+      
+      local_engine.seed( local_seed_type(arma_rng::randi<local_seed_type>()) );
+      
+      for(uword i=0; i<N; ++i)  { mem[i] = eT(local_e_distr(local_engine)); }
+      }
+    #endif
+    }
+  };
+
+
+
 //! @}

include/armadillo_bits/arma_version.hpp

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

include/armadillo_bits/constants.hpp

@@ -26,6 +26,7 @@ namespace priv
     {
     template<typename eT>
     static
+    constexpr
     typename arma_real_only<eT>::result
     nan(typename arma_real_only<eT>::result* junk = nullptr)
       {
@@ -37,6 +38,7 @@ namespace priv
     
     template<typename eT>
     static
+    constexpr
     typename arma_cx_only<eT>::result
     nan(typename arma_cx_only<eT>::result* junk = nullptr)
       {
@@ -50,6 +52,7 @@ namespace priv
     
     template<typename eT>
     static
+    constexpr
     typename arma_integral_only<eT>::result
     nan(typename arma_integral_only<eT>::result* junk = nullptr)
       {
@@ -61,6 +64,7 @@ namespace priv
     
     template<typename eT>
     static
+    constexpr
     typename arma_real_only<eT>::result
     inf(typename arma_real_only<eT>::result* junk = nullptr)
       {
@@ -72,6 +76,7 @@ namespace priv
     
     template<typename eT>
     static
+    constexpr
     typename arma_cx_only<eT>::result
     inf(typename arma_cx_only<eT>::result* junk = nullptr)
       {
@@ -85,6 +90,7 @@ namespace priv
     
     template<typename eT>
     static
+    constexpr
     typename arma_integral_only<eT>::result
     inf(typename arma_integral_only<eT>::result* junk = nullptr)
       {

include/armadillo_bits/diskio_bones.hpp

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

include/armadillo_bits/diskio_meat.hpp

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

include/armadillo_bits/distr_param.hpp

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

include/armadillo_bits/field_meat.hpp

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

include/armadillo_bits/fn_all.hpp

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

include/armadillo_bits/fn_any.hpp

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

include/armadillo_bits/fn_circshift.hpp

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

include/armadillo_bits/fn_cumprod.hpp

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

include/armadillo_bits/fn_cumsum.hpp

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

include/armadillo_bits/fn_det.hpp

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

include/armadillo_bits/fn_diagmat.hpp

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

include/armadillo_bits/fn_diagvec.hpp

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

include/armadillo_bits/fn_diff.hpp

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

include/armadillo_bits/fn_expmat.hpp

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

include/armadillo_bits/fn_flip.hpp

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

include/armadillo_bits/fn_hist.hpp

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

include/armadillo_bits/fn_histc.hpp

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

include/armadillo_bits/fn_index_max.hpp

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

include/armadillo_bits/fn_index_min.hpp

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

include/armadillo_bits/fn_kron.hpp

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

include/armadillo_bits/fn_logmat.hpp

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

include/armadillo_bits/fn_max.hpp

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

include/armadillo_bits/fn_mean.hpp

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

include/armadillo_bits/fn_median.hpp

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

include/armadillo_bits/fn_min.hpp

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

include/armadillo_bits/fn_normalise.hpp

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

include/armadillo_bits/fn_orth_null.hpp

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

include/armadillo_bits/fn_prod.hpp

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

include/armadillo_bits/fn_quantile.hpp

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

include/armadillo_bits/fn_rande.hpp

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

include/armadillo_bits/fn_range.hpp

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

include/armadillo_bits/fn_repelem.hpp

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

include/armadillo_bits/fn_shuffle.hpp

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

include/armadillo_bits/fn_sort.hpp

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

include/armadillo_bits/fn_sort_index.hpp

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

include/armadillo_bits/fn_sqrtmat.hpp

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

include/armadillo_bits/fn_symmat.hpp

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

include/armadillo_bits/glue_hist_meat.hpp

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

include/armadillo_bits/glue_join_bones.hpp

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

include/armadillo_bits/glue_join_meat.hpp

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

include/armadillo_bits/glue_kron_bones.hpp

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

include/armadillo_bits/glue_kron_meat.hpp

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

include/armadillo_bits/glue_polyfit_bones.hpp

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

include/armadillo_bits/glue_polyfit_meat.hpp

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

include/armadillo_bits/glue_polyval_bones.hpp

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

include/armadillo_bits/glue_polyval_meat.hpp

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

include/armadillo_bits/glue_quantile_bones.hpp

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

include/armadillo_bits/glue_quantile_meat.hpp

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

include/armadillo_bits/glue_times_bones.hpp

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

include/armadillo_bits/glue_times_meat.hpp

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

include/armadillo_bits/glue_trapz_bones.hpp

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

include/armadillo_bits/glue_trapz_meat.hpp

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

include/armadillo_bits/mul_gemm.hpp

@@ -59,6 +59,156 @@ struct gemm_emul_tinysq
 
 
 
+struct gemm_emul_large_mp_helper
+  {
+  template<typename eT>
+  arma_hot
+  inline
+  static
+  void
+  copy_row(eT* out_mem, const Mat<eT>& in, const uword row)
+    {
+    const uword n_rows = in.n_rows;
+    const uword n_cols = in.n_cols;
+    
+    const eT* in_mem_row = in.memptr() + row;
+    
+    for(uword i=0; i < n_cols; ++i)
+      {
+      out_mem[i] = (*in_mem_row);
+      
+      in_mem_row += n_rows;
+      }
+    }
+  };
+
+
+
+#if defined(ARMA_USE_OPENMP)
+//! emulation of gemm(), for non-complex matrices only, as it assumes only simple transposes (ie. doesn't do hermitian transposes)
+//! parallelised version
+template<const bool do_trans_A=false, const bool do_trans_B=false, const bool use_alpha=false, const bool use_beta=false>
+struct gemm_emul_large_mp
+  {
+  template<typename eT, typename TA, typename TB>
+  arma_hot
+  inline
+  static
+  void
+  apply
+    (
+          Mat<eT>& C,
+    const TA&      A,
+    const TB&      B,
+    const eT       alpha = eT(1),
+    const eT       beta  = eT(0)
+    )
+    {
+    arma_debug_sigprint();
+    
+    const uword A_n_rows = A.n_rows;
+    const uword A_n_cols = A.n_cols;
+    
+    const uword B_n_rows = B.n_rows;
+    const uword B_n_cols = B.n_cols;
+    
+    if( (do_trans_A == false) && (do_trans_B == false) )
+      {
+      const uword n_threads = uword(mp_thread_limit::get());
+      
+      podarray<eT> tmp(A_n_cols * n_threads, arma_nozeros_indicator());
+      
+      eT* tmp_mem = tmp.memptr();
+      
+      #pragma omp parallel for schedule(static) num_threads(int(n_threads))
+      for(uword row_A=0; row_A < A_n_rows; ++row_A)
+        {
+        const uword thread_id = uword(omp_get_thread_num());
+        
+        eT* A_rowdata = tmp_mem + (A_n_cols * thread_id);
+        
+        gemm_emul_large_mp_helper::copy_row(A_rowdata, A, row_A);
+        
+        for(uword col_B=0; col_B < B_n_cols; ++col_B)
+          {
+          const eT acc = op_dot::direct_dot(B_n_rows, A_rowdata, B.colptr(col_B));
+          
+               if( (use_alpha == false) && (use_beta == false) )  { C.at(row_A,col_B) =       acc;                          }
+          else if( (use_alpha == true ) && (use_beta == false) )  { C.at(row_A,col_B) = alpha*acc;                          }
+          else if( (use_alpha == false) && (use_beta == true ) )  { C.at(row_A,col_B) =       acc + beta*C.at(row_A,col_B); }
+          else if( (use_alpha == true ) && (use_beta == true ) )  { C.at(row_A,col_B) = alpha*acc + beta*C.at(row_A,col_B); }
+          }
+        }
+      }
+    else
+    if( (do_trans_A == true) && (do_trans_B == false) )
+      {
+      const int n_threads = mp_thread_limit::get();
+      
+      #pragma omp parallel for schedule(static) num_threads(n_threads)
+      for(uword col_A=0; col_A < A_n_cols; ++col_A)
+        {
+        // col_A is interpreted as row_A when storing the results in matrix C
+        
+        const eT* A_coldata = A.colptr(col_A);
+        
+        for(uword col_B=0; col_B < B_n_cols; ++col_B)
+          {
+          const eT acc = op_dot::direct_dot(B_n_rows, A_coldata, B.colptr(col_B));
+          
+               if( (use_alpha == false) && (use_beta == false) )  { C.at(col_A,col_B) =       acc;                          }
+          else if( (use_alpha == true ) && (use_beta == false) )  { C.at(col_A,col_B) = alpha*acc;                          }
+          else if( (use_alpha == false) && (use_beta == true ) )  { C.at(col_A,col_B) =       acc + beta*C.at(col_A,col_B); }
+          else if( (use_alpha == true ) && (use_beta == true ) )  { C.at(col_A,col_B) = alpha*acc + beta*C.at(col_A,col_B); }
+          }
+        }
+      }
+    else
+    if( (do_trans_A == false) && (do_trans_B == true) )
+      {
+      Mat<eT> BB;
+      op_strans::apply_mat_noalias(BB, B);
+      
+      gemm_emul_large_mp<false, false, use_alpha, use_beta>::apply(C, A, BB, alpha, beta);
+      }
+    else
+    if( (do_trans_A == true) && (do_trans_B == true) )
+      {
+      // using trans(A)*trans(B) = trans(B*A) equivalency; assuming no hermitian transpose
+      
+      const uword n_threads = uword(mp_thread_limit::get());
+      
+      podarray<eT> tmp(B_n_cols * n_threads, arma_nozeros_indicator());
+      
+      eT* tmp_mem = tmp.memptr();
+      
+      #pragma omp parallel for schedule(static) num_threads(int(n_threads))
+      for(uword row_B=0; row_B < B_n_rows; ++row_B)
+        {
+        const uword thread_id = uword(omp_get_thread_num());
+        
+        eT* B_rowdata = tmp_mem + (B_n_cols * thread_id);
+        
+        gemm_emul_large_mp_helper::copy_row(B_rowdata, B, row_B);
+        
+        for(uword col_A=0; col_A < A_n_cols; ++col_A)
+          {
+          const eT acc = op_dot::direct_dot(A_n_rows, B_rowdata, A.colptr(col_A));
+          
+               if( (use_alpha == false) && (use_beta == false) )  { C.at(col_A,row_B) =       acc;                          }
+          else if( (use_alpha == true ) && (use_beta == false) )  { C.at(col_A,row_B) = alpha*acc;                          }
+          else if( (use_alpha == false) && (use_beta == true ) )  { C.at(col_A,row_B) =       acc + beta*C.at(col_A,row_B); }
+          else if( (use_alpha == true ) && (use_beta == true ) )  { C.at(col_A,row_B) = alpha*acc + beta*C.at(col_A,row_B); }
+          }
+        }
+      }
+    }
+  
+  };
+#endif
+
+
+
 //! emulation of gemm(), for non-complex matrices only, as it assumes only simple transposes (ie. doesn't do hermitian transposes)
 template<const bool do_trans_A=false, const bool do_trans_B=false, const bool use_alpha=false, const bool use_beta=false>
 struct gemm_emul_large
@@ -78,13 +228,28 @@ struct gemm_emul_large
     )
     {
     arma_debug_sigprint();
-
+    
     const uword A_n_rows = A.n_rows;
     const uword A_n_cols = A.n_cols;
     
     const uword B_n_rows = B.n_rows;
     const uword B_n_cols = B.n_cols;
     
+    #if defined(ARMA_USE_OPENMP)
+      {
+      // TODO: replace with more sophisticated threshold mechanism
+      
+      constexpr uword threshold = uword(30);
+      
+      if( (A_n_rows >= threshold) && (A_n_cols >= threshold) && (B_n_rows >= threshold) && (B_n_cols >= threshold) && (mp_thread_limit::in_parallel() == false) )
+        {
+        gemm_emul_large_mp<do_trans_A, do_trans_B, use_alpha, use_beta>::apply(C,A,B,alpha,beta);
+        
+        return;
+        }
+      }
+    #endif
+    
     if( (do_trans_A == false) && (do_trans_B == false) )
       {
       arma_aligned podarray<eT> tmp(A_n_cols);

include/armadillo_bits/mul_gemv.hpp

@@ -203,6 +203,74 @@ struct gemv_emul_helper
 
 
 
+#if defined(ARMA_USE_OPENMP)
+//! Partial emulation of BLAS gemv().
+//! 'y' is assumed to have been set to the correct size (ie. taking into account the transpose)
+//! parallelised version
+template<const bool do_trans_A=false, const bool use_alpha=false, const bool use_beta=false>
+struct gemv_emul_mp
+  {
+  template<typename eT, typename TA>
+  arma_hot
+  inline
+  static
+  void
+  apply( eT* y, const TA& A, const eT* x, const eT alpha = eT(1), const eT beta = eT(0) )
+    {
+    arma_debug_sigprint();
+    
+    const int n_threads = mp_thread_limit::get();
+    
+    const uword A_n_rows = A.n_rows;
+    const uword A_n_cols = A.n_cols;
+    
+    if(do_trans_A == false)
+      {
+      #pragma omp parallel for schedule(static) num_threads(n_threads)
+      for(uword row=0; row < A_n_rows; ++row)
+        {
+        const eT acc = gemv_emul_helper::dot_row_col(A, x, row, A_n_cols);
+        
+             if( (use_alpha == false) && (use_beta == false) )  { y[row] =       acc;               }
+        else if( (use_alpha == true ) && (use_beta == false) )  { y[row] = alpha*acc;               }
+        else if( (use_alpha == false) && (use_beta == true ) )  { y[row] =       acc + beta*y[row]; }
+        else if( (use_alpha == true ) && (use_beta == true ) )  { y[row] = alpha*acc + beta*y[row]; }
+        }
+      }
+    else
+    if(do_trans_A == true)
+      {
+      if(is_cx<eT>::no)
+        {
+        #pragma omp parallel for schedule(static) num_threads(n_threads)
+        for(uword col=0; col < A_n_cols; ++col)
+          {
+          // col is interpreted as row when storing the results in 'y'
+          
+          const eT acc = op_dot::direct_dot(A_n_rows, A.colptr(col), x);
+          
+               if( (use_alpha == false) && (use_beta == false) )  { y[col] =       acc;               }
+          else if( (use_alpha == true ) && (use_beta == false) )  { y[col] = alpha*acc;               }
+          else if( (use_alpha == false) && (use_beta == true ) )  { y[col] =       acc + beta*y[col]; }
+          else if( (use_alpha == true ) && (use_beta == true ) )  { y[col] = alpha*acc + beta*y[col]; }
+          }
+        }
+      else
+        {
+        Mat<eT> AA;
+        
+        op_htrans::apply_mat_noalias(AA, A);
+        
+        gemv_emul_mp<false, use_alpha, use_beta>::apply(y, AA, x, alpha, beta);
+        }
+      }
+    }
+  
+  };
+#endif
+
+
+
 //! \brief
 //! Partial emulation of BLAS gemv().
 //! 'y' is assumed to have been set to the correct size (ie. taking into account the transpose)
@@ -222,6 +290,21 @@ struct gemv_emul
     const uword A_n_rows = A.n_rows;
     const uword A_n_cols = A.n_cols;
     
+    #if defined(ARMA_USE_OPENMP)
+      {
+      // TODO: replace with more sophisticated threshold mechanism
+      
+      constexpr uword threshold = uword(200);
+      
+      if( (A_n_rows >= threshold) && (A_n_cols >= threshold) && (mp_thread_limit::in_parallel() == false) )
+        {
+        gemv_emul_mp<do_trans_A, use_alpha, use_beta>::apply(y, A, x, alpha, beta);
+        
+        return;
+        }
+      }
+    #endif
+    
     if(do_trans_A == false)
       {
       if(A_n_rows == 1)

include/armadillo_bits/op_accu_bones.hpp

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

include/armadillo_bits/op_accu_meat.hpp

@@ -114,7 +114,7 @@ op_accu_mat::apply(const T1& X)
   {
   arma_debug_sigprint();
   
-  if( (is_Mat<T1>::value) || (is_subview_col<T1>::value) || (is_Mat<typename Proxy<T1>::stored_type>::value) || (arma_config::openmp && Proxy<T1>::use_mp) )
+  if( (quasi_unwrap<T1>::has_orig_mem) || (is_Mat<typename Proxy<T1>::stored_type>::value) || (arma_config::openmp && Proxy<T1>::use_mp) )
     {
     const quasi_unwrap<T1> U(X);
     
@@ -602,19 +602,7 @@ op_accu_mat::apply(const subview<eT>& X)
   const uword X_n_rows = X.n_rows;
   const uword X_n_cols = X.n_cols;
   
-  if(X_n_rows == 1)
-    {
-    const Mat<eT>& m = X.m;
-    
-    const uword col_offset = X.aux_col1;
-    const uword row_offset = X.aux_row1;
-    
-    eT val = eT(0);
-    
-    for(uword i=0; i < X_n_cols; ++i)  { val += m.at(row_offset, col_offset + i); }
-    
-    return val;
-    }
+  if(X_n_rows == 1)  { return op_accu_mat::apply( static_cast< const subview_row<eT>& >(X) ); }
   
   if(X_n_cols == 1)  { return arrayops::accumulate( X.colptr(0), X_n_rows ); }
   
@@ -642,6 +630,39 @@ op_accu_mat::apply(const subview_col<eT>& X)
 
 
 
+template<typename eT>
+inline
+eT
+op_accu_mat::apply(const subview_row<eT>& X)
+  {
+  arma_debug_sigprint();  
+  
+  const uword X_m_n_rows = X.m.n_rows;
+  const uword X_n_cols   = X.n_cols;
+  
+  const eT* row_mem = &(X.m.at(X.aux_row1,X.aux_col1));
+  
+  eT val1 = eT(0);
+  eT val2 = eT(0);
+  
+  uword j;
+  
+  for(j=1; j < X_n_cols; j+=2)
+    {
+    val1 += (*row_mem); row_mem += X_m_n_rows;
+    val2 += (*row_mem); row_mem += X_m_n_rows;
+    }
+  
+  if((j-1) < X_n_cols)
+    {
+    val1 += (*row_mem);
+    }
+  
+  return val1 + val2;
+  }
+
+
+
 //
 
 
@@ -926,4 +947,44 @@ op_accu_cube::apply(const CubeToMatOp<T1, op_omit_cube>& in)
 
 
 
+template<typename eT>
+inline
+eT
+op_accu_cube::apply(const subview_cube<eT>& sv)
+  {
+  arma_debug_sigprint();  
+  
+  if(sv.n_elem == 0)  { return eT(0); }
+  
+  const uword sv_nr = sv.n_rows;
+  const uword sv_nc = sv.n_cols;
+  const uword sv_ns = sv.n_slices;
+  
+  eT acc = eT(0);
+  
+  if( (sv_nr == 1)  && (sv_nc == 1) && (sv.aux_slice1 == 0) )
+    {
+    const uword sv_m_n_elem_slice = sv.m.n_elem_slice;
+    
+    const eT* sv_m_ptr = &( sv.m.at(sv.aux_row1, sv.aux_col1, 0) );
+    
+    for(uword s=0; s < sv_ns; ++s)
+      {
+      acc += (*sv_m_ptr);  sv_m_ptr += sv_m_n_elem_slice;
+      }
+    }
+  else
+    {
+    for(uword s=0; s < sv_ns; ++s)
+    for(uword c=0; c < sv_nc; ++c)
+      {
+      acc += arrayops::accumulate(sv.slice_colptr(s,c), sv_nr);
+      }
+    }
+  
+  return acc;
+  }
+
+
+
 //! @}

include/armadillo_bits/op_all_meat.hpp

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

include/armadillo_bits/op_any_meat.hpp

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

include/armadillo_bits/op_circshift_bones.hpp

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

include/armadillo_bits/op_circshift_meat.hpp

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

include/armadillo_bits/op_clamp_bones.hpp

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

include/armadillo_bits/op_clamp_meat.hpp

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

include/armadillo_bits/op_cumprod_bones.hpp

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

include/armadillo_bits/op_cumprod_meat.hpp

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

include/armadillo_bits/op_cumsum_bones.hpp

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

include/armadillo_bits/op_cumsum_meat.hpp

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

include/armadillo_bits/op_diff_bones.hpp

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

include/armadillo_bits/op_diff_meat.hpp

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

include/armadillo_bits/op_flip_bones.hpp

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

include/armadillo_bits/op_flip_meat.hpp

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

include/armadillo_bits/op_htrans_bones.hpp

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

include/armadillo_bits/op_htrans_meat.hpp

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

include/armadillo_bits/op_max_bones.hpp

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

include/armadillo_bits/op_max_meat.hpp

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

include/armadillo_bits/op_mean_bones.hpp

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

include/armadillo_bits/op_mean_meat.hpp

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

include/armadillo_bits/op_median_bones.hpp

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

include/armadillo_bits/op_median_meat.hpp

@@ -53,6 +53,25 @@ op_median::apply(Mat<typename T1::elem_type>& out, const Op<T1,op_median>& expr)
 
 
 
+template<typename T1>
+inline
+void
+op_median::apply(Mat_noalias<typename T1::elem_type>& out, const Op<T1,op_median>& expr)
+  {
+  arma_debug_sigprint();
+  
+  const quasi_unwrap<T1> U(expr.m);
+  
+  const uword dim = expr.aux_uword_a;
+  
+  arma_conform_check( U.M.internal_has_nan(), "median(): detected NaN"                   );
+  arma_conform_check( (dim > 1),              "median(): parameter 'dim' must be 0 or 1" );
+  
+  op_median::apply_noalias(out, U.M, dim);
+  }
+
+
+
 template<typename eT>
 inline
 void

include/armadillo_bits/op_min_bones.hpp

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

include/armadillo_bits/op_min_meat.hpp

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

include/armadillo_bits/op_nonzeros_bones.hpp

@@ -32,6 +32,9 @@ struct op_nonzeros
   
   template<typename T1>
   static inline void apply(Mat<typename T1::elem_type>& out, const Op<T1, op_nonzeros>& X);
+  
+  template<typename T1>
+  static inline void apply(Mat_noalias<typename T1::elem_type>& out, const Op<T1, op_nonzeros>& X);
   };
 
 

include/armadillo_bits/op_nonzeros_meat.hpp

@@ -100,4 +100,20 @@ op_nonzeros::apply(Mat<typename T1::elem_type>& out, const Op<T1, op_nonzeros>&
 
 
 
+template<typename T1>
+inline
+void
+op_nonzeros::apply(Mat_noalias<typename T1::elem_type>& out, const Op<T1, op_nonzeros>& X)
+  {
+  arma_debug_sigprint();
+  
+  const Proxy<T1> P(X.m);
+  
+  if(P.get_n_elem() == 0)  { out.set_size(0,1); return; }
+  
+  op_nonzeros::apply_noalias(out, P);
+  }
+
+
+
 //! @}

include/armadillo_bits/op_repmat_bones.hpp

@@ -28,6 +28,8 @@ struct op_repmat
   template<typename obj> inline static void apply_noalias(Mat<typename obj::elem_type>& out, const obj& X, const uword copies_per_row, const uword copies_per_col);
   
   template<typename T1> inline static void apply(Mat<typename T1::elem_type>& out, const Op<T1,op_repmat>& in);
+  
+  template<typename T1> inline static void apply(Mat_noalias<typename T1::elem_type>& out, const Op<T1,op_repmat>& in);
   };