---
title: "Functions of vectors, matrices, and cubes"
output: rmarkdown::html_vignette
bibliography: "references.bib"
vignette: >
  %\VignetteIndexEntry{Functions of vectors, matrices, and cubes}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
```

**This vignette is adapted from the official Armadillo
[documentation](https://arma.sourceforge.net/docs.html).**

# Contents

| Function                | Description                                                                 |
|-------------------------|-----------------------------------------------------------------------------|
| `abs`                   | [Obtain magnitude of each element](#abs)                                    |
| `accu`                  | [Accumulate (sum) all elements](#accu)                                      |
| `affmul`                | [Affine matrix multiplication](#affmul)                                     |
| `all`                   | [Check whether all elements are non-zero, or satisfy a relational condition](#all) |
| `any`                   | [Check whether any element is non-zero, or satisfies a relational condition](#any) |
| `approx_equal`          | [Approximate equality](#approx_equal)                                       |
| `arg`                   | [Phase angle of each element](#arg)                                         |
| `as_scalar`             | [Convert 1x1 matrix to pure scalar](#as_scalar)                             |
| `clamp`                 | [Obtain clamped elements according to given limits](#clamp)                 |
| `cond`                  | [Condition number of matrix](#cond)                                         |
| `conj`                  | [Obtain complex conjugate of each element](#conj)                           |
| `conv_to`               | [Convert/cast between matrix types](#conv_to)                               |
| `cross`                 | [Cross product](#cross)                                                     |
| `cumsum`                | [Cumulative sum](#cumsum)                                                   |
| `cumprod`               | [Cumulative product](#cumprod)                                              |
| `det`                   | [Determinant](#det)                                                         |
| `diagmat`               | [Generate diagonal matrix from given matrix or vector](#diagmat)            |
| `diagvec`               | [Extract specified diagonal](#diagvec)                                      |
| `diags`                 | [Generate a dense matrix with diagonals specified by column vectors](#diags) |
| `diff`                  | [Differences between adjacent elements](#diff)                              |
| `dot/cdot/norm_dot`     | [Dot product](#dot)                                         |
| `eps`                   | [Obtain distance of each element to next largest floating point representation](#eps) |
| `expmat`                | [Matrix exponential](#expmat)                                               |
| `expmat_sym`            | [Matrix exponential of symmetric matrix](#expmat_sym)                       |
| `find`                  | [Find indices of non-zero elements, or elements satisfying a relational condition](#find) |
| `find_finite`           | [Find indices of finite elements](#find_finite)                             |
| `find_nonfinite`        | [Find indices of non-finite elements](#find_nonfinite)                      |
| `find_nan`              | [Find indices of NaN elements](#find_nan)                                   |
| `find_unique`           | [Find indices of unique elements](#find_unique)                             |
| `fliplr/flipud`         | [Flip matrix left to right or upside down](#fliplr)                         |
| `imag/real`             | [Extract imaginary/real part](#imag)                                        |
| `ind2sub`               | [Convert linear index to subscripts](#ind2sub)                              |
| `index_min / index_max` | [Indices of extremum values](#index_min)                                    |
| `inplace_trans`         | [In-place dense transpose](#inplace_trans)                                  |
| `intersect`             | [Find common elements in two vectors/matrices](#intersect)                  |
| `join_rows / join_cols` | [Concatenation of matrices](#join_rows)                                     |
| `join_slices`           | [Concatenation of cubes](#join_slices)                                      |
| `kron`                  | [Kronecker tensor product](#kron)                                           |
| `log_det`               | [Log determinant](#log_det)                                                 |
| `log_det_sympd`         | [Log determinant of symmetric positive definite matrix](#log_det_sympd)     |
| `logmat`                | [Matrix logarithm](#logmat)                                                 |
| `logmat_sympd`          | [Matrix logarithm of symmetric matrix](#logmat_sympd)                       |
| `min/max`               | [Return extremum values](#min)                                              |
| `nonzeros`              | [Return non-zero values](#nonzeros)                                         |
| `norm`                  | [Various norms of vectors and matrices](#norm)                              |
| `norm2est`              | [Fast estimate of the matrix 2-norm](#norm2est)                             |
| `normalise`             | [Normalise vectors to unit p-norm](#normalise)                              |
| `pow`                   | [Element-wise power](#pow)                                                  |
| `powmat`                | [Matrix power](#powmat)                                                     |
| `prod`                  | [Product of elements](#prod)                                                |
| `rank`                  | [Rank of matrix](#rank)                                                     |
| `rcond`                 | [Reciprocal condition number](#rcond)                                       |
| `repelem`               | [Replicate elements](#repelem)                                              |
| `repmat`                | [Replicate matrix in block-like fashion](#repmat)                           |
| `reshape`               | [Change size while keeping elements](#reshape)                              |
| `resize`                | [Change size while keeping elements and preserving layout](#resize)         |
| `reverse`               | [Reverse order of elements](#reverse)                                       |
| `roots`                 | [Roots of polynomial](#roots)                                               |
| `shift`                 | [Shift elements](#shift)                                                    |
| `shuffle`               | [Randomly shuffle elements](#shuffle)                                       |
| `size`                  | [Obtain dimensions of given object](#size)                                  |
| `sort`                  | [Sort elements](#sort)                                                      |
| `sort_index`            | [Vector describing sorted order of elements](#sort_index)                   |
| `spdiags`               | [Generate a sparse matrix with diagonals specified by column vectors](#spdiags) |
| `sqrtmat`               | [Square root of matrix](#sqrtmat)                                           |
| `sqrtmat_sympd`         | [Square root of symmetric matrix](#sqrtmat_sympd)                           |
| `sum`                   | [Sum of elements](#sum)                                                     |
| `sub2ind`               | [Convert subscripts to linear index](#sub2ind)                              |
| `symmatu / symmatl`     | [Generate symmetric matrix from given matrix](#symmatu-symmatl)            |
| `trace`                 | [Sum of diagonal elements](#trace)                                          |
| `trans / strans`        | [Transpose of matrix](#trans-strans)                                               |
| `trapz`                 | [Trapezoidal numerical integration](#trapz)                                 |
| `trimatu / trimatl`     | [Copy upper/lower triangular part](#trimatu-trimatl)                       |
| `trimatu_ind / trimatl_ind` | [Obtain indices of upper/lower triangular part](#trimatu_ind-trimatl_ind) |
| `unique`                | [Return unique elements](#unique)                                           |
| `vecnorm`               | [Obtain vector norm of each row or column of a matrix](#vecnorm)            |
| `vectorise`             | [Flatten matrix into vector](#vectorise)                                    |
| `misc functions`        | [Miscellaneous element-wise functions: exp, log, sqrt, round, sign, ...](#misc-functions) |
| `trig functions`        | [Trigonometric element-wise functions: cos, sin, tan, ...](#trig-functions) |

# Absolute value {#abs}

The `abs()` function computes the absolute value of each element in a vector,
matrix, or cube.

Usage:

```cpp
Y = abs(X); // for non-complex X
real_object_type Y = abs(X); // for complex X
```

For the non-complex case, `X` and `Y` must have the same type, such as mat or
cube.

For the complex case, `Y` must be the real counterpart to the type of `X`. If
`X` has the type `cx_mat`, then the type of `Y` must be `mat `.

## Examples

```cpp
[[cpp11::register]] doubles_matrix<> abs1_(const int& n) {
  mat A(n, n, fill::randu);
  mat B = abs(A);

  cx_mat X(n, n, fill::randu);
  mat Y = abs(X);

  mat res = B + Y;

  return as_doubles_matrix(res);
}
```

#  Accumulate (sum) all elements {#accu}

The `accu()` function computes the sum of all elements in a vector, matrix, or
cube.

## Examples

```cpp
[[cpp11::register]] double accu1_(const int& n) {
  mat A(n, n, fill::randu);
  mat B(n, n, fill::randu);

  double x = accu(A);

  // accu(A % B) is a "multiply-and-accumulate" operation
  // as operator % performs element-wise multiplication
  double y = accu(A % B);

  return (x + y);
}
```

# Affine matrix multiplication {#affmul}

The `affmul()` function computes matrix multiplication for `A` and `B` with an
extended form of `B`. `A` is typically an affine transformation matrix. `B` can
be a vector or matrix, and is treated as having an additional row of ones.

The number of columns in `A` must be equal to number of rows in the extended
form of `B` (e.g., `A.n_cols = B.n_rows + 1`).

If `A`has dimensions 3x3 and `B` 2x1, the equivalent matrix multiplication is:

```
⎡ C0 ⎤   ⎡ A00 A01 A02 ⎤   ⎡ B0 ⎤
⎢ C1 ⎥ = ⎢ A10 A11 A12 ⎥ x ⎢ B1 ⎥
⎣ C2 ⎦   ⎣ A20 A21 A22 ⎦   ⎣ 1  ⎦
```

If `A` has dimensions 2x3 and `B` 2x1, the equivalent matrix multiplication is:

```
⎡ C0 ⎤   ⎡ A00 A01 A02 ⎤   ⎡ B0 ⎤
⎢ C1 ⎥ = ⎢ A10 A11 A12 ⎥ x ⎢ B1 ⎥
                           ⎣ 1  ⎦
```

## Examples

```cpp
[[cpp11::register]] doubles affmul1_(const int& n) {
  mat A(n, n + 1, fill::randu);
  vec B(n, fill::randu);

  vec C = affmul(A, B);

  return as_doubles(C);
}
```

# Check whether all elements are non-zero, or satisfy a relational condition {#all}

The `all()` function checks whether all elements in a vector, matrix or cube
are non-zero, or satisfy a relational condition. It returns true/false booleans
for vectors and 0/1 vectors for matrices to indicate if the condition is met
for each row or column.

Usage:

```cpp
all(vector);
all(matrix);
all(matrix, dimension); // dimension = 0 -> returns a row vector urowvec/umat
                        // dimension = 1 -> returns a column vector ucolvec/umat
```

## Examples

```cpp
[[cpp11::register]] logicals all1_(const int& n) {
  vec V(n, fill::randu);
  mat X(n, n, fill::randu);

  // true if vector V has all non-zero elements
  bool status1 = all(V);

  // true if vector V has all elements greater than 0.5
  bool status2 = all(V > 0.5);

  // true if matrix X has all elements greater than 0.6;
  // note the use of vectorise()
  bool status3 = all(vectorise(X) > 0.6);

  // row vector indicating which columns of X have all elements greater than 0.7
  umat A = all(X > 0.7);

  writable::logicals res(4);
  res[0] = status1;
  res[1] = status2;
  res[2] = status3;
  res[3] = all(vectorise(A) == 1);  // true if all elements of A are 1

  return res;
}
```

# Check whether any element is non-zero, or satisfies a relational condition {#any}

The `any()` function checks whether any element in a vector, matrix or cube is
non-zero, or satisfies a relational condition. It returns true/false booleans
for vectors and 0/1 vectors for matrices to indicate if the condition is met
for any row or column.

Usage:

```cpp
any(vector);
any(matrix);
any(matrix, dimension); // dimension = 0 -> returns a row vector urowvec/umat
                        // dimension = 1 -> returns a column vector ucolvec/umat
```

## Examples

```cpp
[[cpp11::register]] logicals any1_(const int& n) {
  vec V(n, fill::randu);
  mat X(n, n, fill::randu);

  // true if vector V has any non-zero elements
  bool status1 = any(V);

  // true if vector V has any elements greater than 0.5
  bool status2 = any(V > 0.5);

  // true if matrix X has any elements greater than 0.6;
  // note the use of vectorise()
  bool status3 = any(vectorise(X) > 0.6);

  // row vector indicating which columns of X have any elements greater than 0.7
  umat A = any(X > 0.7);

  writable::logicals res(4);
  res[0] = status1;
  res[1] = status2;
  res[2] = status3;
  res[3] = any(vectorise(A) == 1);  // true if any element of A is 1

  return res;
}
```

# Approximate equality {#approx_equal}

The `approx_equal()` function checks whether two vectors, matrices or cubes are
approximately equal. It returns true if all corresponding elements have
differences less than or equal to a given tolerance.

Usage:

```cpp
approx_equal(A, B, method, tol)
approx_equal(A, B, method, abs_tol, rel_tol)
```

The `method` parameter specifies the method used to compare the elements:

* `method = "absdiff"`: absolute difference (e.g., `|A - B| <= tol`)
* `method = "reldiff"`: relative difference (e.g., `|A - B| / max(|A|, |B|) <= tol`)
* `method = "both"`: absolute or relative difference (e.g., `|A - B| <= tol || |A - B| / max(|A|, |B|) <= tol`)

## Examples

```cpp
[[cpp11::register]] bool approx_equal1_(const int& n) {
  mat A(n, n, fill::randu);
  mat B = A + 0.001;

  bool same1 = approx_equal(A, B, "absdiff", 0.002);

  mat C = 1000 * randu<mat>(n, n);
  mat D = C + 1;

  bool same2 = approx_equal(C, D, "reldiff", 0.1);

  bool same3 = approx_equal(C, D, "both", 2, 0.1);

  bool all_same = same1 && same2 && same3;

  return all_same;
}
```

# Phase angle of each element {#arg}

The `arg()` function computes the phase angle of each element in a vector,
matrix or cube. For non-complex elements, the input is treated as a complex
element with zero imaginary component. For complex elements, the input must be
of the same and the output the real counterpart type.

Usage:

```cpp
real_object_type Y = arg(X);
```

## Examples

```cpp
[[cpp11::register]] doubles_matrix<> arg1_(const int& n) {
  cx_mat X(n, n, fill::randu);
  mat Y = arg(X);

  return as_doubles_matrix(Y);
}
```

# Convert 1x1 matrix to pure scalar {#as_scalar}

The `as_scalar()` function converts a 1x1 matrix to a scalar (e.g.,
`double/int`). It is useful when you want to extract a single element from a
matrix or an operation (e.g., converting the result of a dot/inner product to a
scalar).

## Examples

```cpp
[[cpp11::register]] double as_scalar1_(const int& n) {
  rowvec r(n, fill::randu);
  colvec q(n, fill::randu);

  mat X(n, n, fill::randu);

  // examples of expressions which have optimised implementations
  double a = as_scalar(r*q);
  double b = as_scalar(r*X*q);
  double c = as_scalar(r*diagmat(X)*q);
  double d = as_scalar(r*inv(diagmat(X))*q);

  return (a + b + c + d);
}
```

# Obtain clamped elements according to given limits {#clamp}

The `clamp()` function clamps each element in a vector, matrix or cube to a
given range. Any value less than the lower limit is set to the lower limit, and
any value greater than the upper limit is set to the upper limit.

For objects with complex elements, the real and imaginary components are
clamped separately.

If the input is a sparse matrix, only the non-zero elements are clamped.

## Example

```cpp
[[cpp11::register]] doubles_matrix<> clamp1_(const int& n) {
  mat A(n, n, fill::randu);
  mat B = clamp(A, 0.2, 0.8);
  mat C = clamp(A, A.min(), 0.8);
  mat D = clamp(A, 0.2, A.max());

  mat res = B + C + D;

  return as_doubles_matrix(res);
}
```

# Condition number of matrix {#cond}

The `cond()` function computes the condition number of a matrix. The condition
number is the ratio of the largest singular value to the smallest singular
value. It is a measure of how well the matrix can be inverted, a matrix with a
value close to 1 is well-conditioned, and a matrix with a large value is
ill-conditioned. The computation is based on the singular value decomposition.

## Examples

```cpp
[[cpp11::register]] double cond1_(const int& n) {
  mat A(n, n);
  A.eye(); // the identity matrix has a condition number of 1

  double cond_num = cond(A);

  return cond_num;
}
```

## Caveat

Calculating the approximate reciprocal condition number via `rcond()` is
considerably more efficient.

# Obtain complex conjugate of each element {#conj}

The `conj()` function computes the complex conjugate of each element in a
complex matrix or cube.

## Examples

```cpp
[[cpp11::register]] list conj1_(const int& n) {
  cx_mat X(n, n, fill::randu);
  cx_mat Y = conj(X);
  return as_complex_matrix(Y);
}
```

# Convert/cast between matrix types {#conv_to}

The `conv_to()` function converts a matrix or cube to a different type. It can
convert `mat` to `imat`, `cube` to `icube`, `mat` into `colvec` or any other
casting that preserves data (e.g., a matrix that cannot be interpreted as a
vector is not a valid casting). It can also be used to convert a matrix/vector
into a `std::vector` object.

Usage:

```cpp
conv_to<type>::from(X) 
```

## Examples

```cpp
[[cpp11::register]] doubles conv_to1_(const int& n) {
  mat A(n, n, fill::randu);
  fmat B = conv_to<fmat>::from(A);

  std::vector<double> x(B.n_elem);

  int i, N = static_cast<int>(B.n_elem);
  for (i = 0; i < N; ++i) { x[i] = B(i); }

  colvec y = conv_to<colvec>::from(x);
  std::vector<double> z = conv_to<std::vector<double>>::from(y);

  return as_doubles(z);
}
```

## Caveat

To convert an expression that results in a 1x1 matrix to a pure scalar value,
use `as_scalar()`.

# Cross product {#cross}

The `cross()` function computes the cross product of two vectors under the
assumption that the vectors are three-dimensional.

## Examples

```cpp
[[cpp11::register]] doubles cross1_(const int& n) {
  vec A(n, fill::randu);
  vec B(n, fill::randu);

  vec C = cross(A, B);

  return as_doubles(C);
}
```

# Cumulative sum {#cumsum}

The `cumsum()` function computes the cumulative sum of elements in a vector or
matrix. For a vector, it returns a vector of the same orientation. For a matrix,
it returns a matrix with the cumulative sum along the specified dimension (the
default is along columns with `dimension = 0`).

Usage:

```cpp
cumsum(vector);
cumsum(matrix, dimension); // dimension = 0 -> cumulative sum along columns
                           // dimension = 1 -> cumulative sum along rows
```

## Examples

```cpp
[[cpp11::register]] doubles cumsum1_(const int& n) {
  mat A(n, n, fill::randu);
  mat B = cumsum(A);
  mat C = cumsum(A, 1);

  vec x(n, fill::randu);
  vec y = cumsum(x);

  writable::doubles res(3);
  res[0] = accu(B);
  res[1] = accu(C);
  res[2] = accu(y);
  
  return res;
}
```

# Cumulative product {#cumprod}

The `cumprod()` function computes the cumulative product of elements in a vector
or matrix. For a vector, it returns a vector of the same orientation. For a
matrix, it returns a matrix with the cumulative product along the specified
dimension (the default is along columns with `dimension = 0`).

Usage:

```cpp
cumprod(vector);
cumprod(matrix, dimension); // dimension = 0 -> cumulative prod along columns
                            // dimension = 1 -> cumulative prod along rows
```

## Examples

```cpp
[[cpp11::register]] doubles cumprod1_(const int& n) {
  mat A(n, n, fill::randu);
  mat B = cumprod(A);
  mat C = cumprod(A, 1);

  vec x(n, fill::randu);
  vec y = cumprod(x);

  writable::doubles res(3);
  res[0] = accu(B);
  res[1] = accu(C);
  res[2] = accu(y);
  
  return res;
}
```

# Determinant {#det}

The `det()` function computes the determinant of a square matrix. It is based
on the LU decomposition. If the input is a not a square matrix, the function
throws a `std::runtime_error` exception.

Usage:

```
val = det(X); // store a scalar
det(val, A); // store the determinant in val and return true if successful
```

If the calculation fails:

* `val = det(A)` throws a `std::runtime_error` exception
* `det(val,A)` returns a bool set to false (exception is not thrown)

## Examples

```cpp
[[cpp11::register]] doubles det1_(const int& n) {
  mat A(n, n, fill::randu);
  double val1 = det(A);

  double val2;
  mat B(n, n, fill::randu);
  bool success2 = det(val2, B);

  return writable::doubles({val1, val2, static_cast<double>(success2)});
}
```

# Generate diagonal matrix from given matrix or vector {#diagmat}

The `diagmat()` function generates a diagonal matrix from a given vector or
matrix. If the input is a vector, the output is a square matrix with the vector
as the diagonal. If the input is a matrix, the output is a square matrix with
the diagonal elements from the input matrix. Any element outside the diagonal
is set to zero. The default is the main diagonal (`k = 0`).

Usage:

```cpp
diagmat(vector);
diagmat(matrix);
diagmat(matrix, k); // k = 0 -> main diagonal
                    // k > 0 -> above main diagonal
                    // k < 0 -> below main diagonal
```

## Examples

```cpp
[[cpp11::register]] doubles_matrix<> diagmat1_(const int& n) {
  mat A(n, n, fill::randu);
  mat B = diagmat(A);
  mat C = diagmat(A, 1);

  vec v(n, fill::randu);
  mat D = diagmat(v); // NxN diagonal matrix
  mat E = diagmat(v, 1); // (N+1)x(N+1) diagonal matrix

  mat res = B + C + D;  
  res += E.submat(0, 0, 1, 1); // the result is an upper triangular matrix

  return as_doubles_matrix(res);
}
```

# Extract specified diagonal {#diagvec}

The `diagvec()` function extracts the specified diagonal from a matrix. The
default is the main diagonal (`k = 0`).

Usage:

```cpp
diagvec(matrix);
diagvec(matrix, k); // k = 0 -> main diagonal
                    // k > 0 -> above main diagonal
                    // k < 0 -> below main diagonal
```

## Examples

```cpp
[[cpp11::register]] doubles diagvec1_(const int& n) {
  mat A(n, n, fill::randu);
  vec B = diagvec(A);
  vec C = diagvec(A, 1);

  vec res = B.subvec(0, 1) + C;

  return as_doubles(res);
}
```

# Generate a dense matrix with diagonals specified by column vectors {#diags}

The `diags()` function generates a dense matrix with diagonals specified by
column vectors from an input matrix and a vector to indicate the diagonals.

Usage:

```cpp
diags(matrix, vector, number_of_rows, number_of_columns);
```

Each element in the input vector specifies diagonal `k`, where:

* `k = 0` is the main diagonal
* `k > 0` is above the main diagonal
* `k < 0` is below the main diagonal

## Examples

```cpp
[[cpp11::register]] doubles_matrix<> diags1_(const int& n) {
  mat V(n, n, fill::randu);
  ivec D = {0, -1};
  mat X = diags(V, D, n, n); // lower triangular matrix
  return as_doubles_matrix(X);
}
```

# Differences between adjacent elements {#diff}

The `diff()` function computes the differences between adjacent elements in a
vector or matrix. For a vector, the output is a vector of length `n-k` (the
default is `k = 1`). For a matrix, the output is a matrix with `n-k` rows
when `dim = 0` (the default) and `m-k` columns when `dim = 1`. If `k` is greater
than the length of the vector or the number or rows/columns, the output is an
empty vector/matrix.

Usage:

```cpp
diff(vector);
diff(vector, k);

diff(matrix);
diff(matrix, k);
diff(matrix, k, dim); // dim = 0 -> differences along columns
                      // dim = 1 -> differences along rows
```

## Examples

```cpp
[[cpp11::register]] doubles_matrix<> diff1_(const int& n) {
  vec a = randu<vec>(n);
  vec b = diff(a);

  mat res(n, 2, fill::zeros);
  
  res.col(0) = a;

  for (int i = 1; i < n; ++i) {
    res(i, 1) = b(i - 1);
  }

  return as_doubles_matrix(res);
}
```

# Dot product {#dot}

The `dot()`, `cdot()`, and `norm_dot()` functions compute the dot product of two
vectors. The `cdot()` function computes the complex conjugate dot product, and
the `norm_dot()` function computes the dot product and normalises the result by
the product of the Euclidean norms of the input vectors.

## Examples

```cpp
[[cpp11::register]] doubles dot1_(const int& n) {
  vec A(n, fill::randu);
  vec B(n, fill::randu);
  return writable::doubles({dot(A, B), cdot(A, B), norm_dot(A, B)});
}
```

## Caveat

`norm()` is more robust for calculating the norm, as it handles underflows and
overflows.

# Obtain distance of each element to next largest floating point representation {#eps}

The `eps()` function computes the distance of each element in a scalar, vector
or  matrix to the next largest floating point representation. For vector input,
the output is a vector of the same orientation and length. For matrix input,
the output is a matrix of the same dimensions.

## Examples

```cpp
[[cpp11::register]] doubles_matrix<> eps1_(const int& n) {
  mat A(n, n, fill::randu);
  mat B = eps(A);
  return as_doubles_matrix(B);
}
```

# Matrix exponential {#expmat}

The `expmat()` function computes the matrix exponential of a square matrix. If
the matrix exponential cannot be computed, the function throws a
`std::runtime_error`, same if the input is not a square matrix.

## Examples

```cpp
[[cpp11::register]] doubles_matrix<> expmat1_(const int& n) {
  mat A(n, n, fill::randu);
  mat B = expmat(A);
  return as_doubles_matrix(B);
}
```

## Caveats

* The matrix exponential operation is generally not the same as applying the
  `exp()` function to each element.
* If the input matrix is symmetric, `expmat_sym()` is faster.

# Matrix exponential of symmetric matrix {#expmat_sym}

The `expmat_sym()` function computes the matrix exponential of a symmetric or
Hermitian matrix. If the matrix exponential cannot be computed, the function
throws a `std::runtime_error`, same if the input is not a square matrix.

## Examples

```cpp
[[cpp11::register]] doubles_matrix<> expmat_sym1_(const int& n) {
  mat A(n, n, fill::randu);
  A = A + A.t(); // make A symmetric
  mat B = expmat_sym(A);
  return as_doubles_matrix(B);
}
```

# Find indices of non-zero elements, or elements satisfying a relational condition {#find}

The `find()` function returns the indices of non-zero elements in a vector,
or that satisfy a relational condition in a vector or matrix. The output is a
vector of indices (`uvec`).

Usage:

```cpp
find(vector);
find(vector, k);
find(vector, k, s);

find(matrix);
find(matrix, k);
find(matrix, k, s);
```

The parameter `k` (`k=0` by default) returns the indices of all non-zero
elements or elements that meet the condition. The optional parameter
`s = "first"` returns the first `m` non-zero indices or indices that meet the
condition, and `s = "last"` returns the last `m` non-zero indices or indices
that meet the condition.

## Examples

```cpp
[[cpp11::register]] list find1_(const int& n) {
  mat A(n, n, fill::randu);
  mat B(n, n, fill::randu);

  uvec q1 = find(A > B);
  uvec q2 = find(A > 0.5);
  uvec q3 = find(A > 0.5, 3, "last");

  // change elements of A greater than 0.5 to 1
  A.elem(find(A > 0.5)).ones();

  return writable::list(as_integers(q1), as_integers(q2), as_integers(q3));
}
```

## Caveats

* To clamp values to an interval, `clamp()` is more efficient.
* To replace a specific value, `.replace()` is more efficient.

# Find indices of finite elements {#find_finite}

The `find_finite()` function returns the indices of finite elements in a vector
or matrix. The output is a vector of indices (`uvec`).

## Examples

```cpp
[[cpp11::register]] integers find_finite1_(const int& n) {
  mat A(n, n, fill::randu);
  uvec q = find_finite(A);
  return as_integers(q);
}
```

# Find indices of non-finite elements {#find_nonfinite}

The `find_nonfinite()` function returns the indices of non-finite elements in a
vector or matrix. The output is a vector of indices (`uvec`).

## Examples

```cpp
[[cpp11::register]] integers find_nonfinite1_(const int& n) {
  mat A(n, n, fill::randu);
  A(0, 0) = datum::inf;
  uvec q = find_nonfinite(A);
  return as_integers(q);
}
```

## Caveat

To replace instances of a specific non-finite value (eg. `NaN` or `Inf`), it is
more efficient to use `.replace()`.

# Find indices of NaN elements {#find_nan}

The `find_nan()` function returns the indices of NaN elements in a vector or
matrix. The output is a vector of indices (`uvec`).

## Examples

```cpp
[[cpp11::register]] integers find_nan1_(const int& n) {
  mat A(n, n, fill::randu);
  A(0, 0) = datum::nan;
  uvec q = find_nan(A);
  return as_integers(q);
}
```

## Caveat

To replace instances of `NaN` values, it is more efficient to use `.replace()`.

# Find indices of unique elements {#find_unique}

The `find_unique()` function returns the indices of unique elements in a vector
or matrix. The output is a vector of indices (`uvec`).

## Examples

```cpp
[[cpp11::register]] integers find_unique1_(const int& n) {
  mat A(n, n, fill::randu);
  A(0, 0) = A(1, 1);
  uvec q = find_unique(A);
  return as_integers(q);
}
```

# Flip matrix left to right or upside down {#fliplr}

The `fliplr()` function generates a copy of the input matrix with the order of
the columns reversed, and the `flipud()` function generates a copy of the input
matrix with the order of the rows reversed.

## Examples

```cpp
[[cpp11::register]] list flip1_(const int& n) {
  mat A(n, n, fill::randu);
  mat B = fliplr(A);
  mat C = flipud(A);

  writable::list res(3);
  res[0] = as_doubles_matrix(A);
  res[1] = as_doubles_matrix(B);
  res[2] = as_doubles_matrix(C);

  return res;
}
```

# Extract imaginary/real part {#imag}

The `imag()` and `real()` functions extract the imaginary and real parts of
each element in a complex matrix, respectively.

## Examples

```cpp
[[cpp11::register]] list imag1_(const int& n) {
  cx_mat X(n, n, fill::randu);
  mat Y = imag(X);
  mat Z = real(X);

  writable::list res(2);
  res[0] = as_doubles_matrix(Y);
  res[1] = as_doubles_matrix(Z);

  return res;
}
```

## Caveat

To convert a complex matrix to a list of real matrices, it is more efficient to
use `as_complex_matrix()`.

# Convert linear index to subscripts {#ind2sub}

The `ind2sub()` function converts a linear index or vector of indexes to
subscripts. The output is a vector of indices (`uvec`) if the input index is a
scalar, and a matrix of indices (`umat`) if the input index is a vector.

Usage:

```cpp
uvec sub = ind2sub(size(X), index)
uvec sub = ind2sub(size(n_rows, n_cols), index)
uvec sub = ind2sub(size(n_rows, n_cols, n_slices), index)

umat sub = ind2sub(size(X), vector_of_indices)
umat sub = ind2sub(size(n_rows, n_cols), vector_of_indices)
umat sub = ind2sub(size(n_rows, n_cols, n_slices), vector_of_indices)
```

## Examples

```cpp
[[cpp11::register]] list ind2sub1_(const int& n) {
  mat M(n, n, fill::randu);

  uvec s = ind2sub(size(M), n);

  uvec indices = find(M > 0.5);
  umat t       = ind2sub(size(M), indices);

  cube Q(2,3,4);

  uvec u = ind2sub(size(Q), 8);

  writable::list res(3);
  res[0] = as_integers(s);
  res[1] = as_integers_matrix(t);
  res[2] = as_integers(u);

  return res;
}
```

# Indices of extremum values {#index_min}

The `index_min()` and `index_max()` functions return the indices of the minimum
and maximum values in a vector, matrix or cube. For an input vector, the output
is a scalar index (`uword`). For an input matrix, the output is a vector of
indices (`uvec`) with row orientation for the argument `dim = 0` (default) with
the min/max for each column, and column orientation for `dim = 1` with the
min/max for each row. For an input cube, the output is a cube of indices
(`ucube`) with the min/max for each slice's columns when `dim = 0`, the min/max
for each slice's rows when `dim = 1`, and the min/max for each slice when
`dim = 2`. For complex objects, the absolute value is used to compare the
elements.

Usage:

```cpp
// index_max is analogous

index_min(vector)

index_min(matrix)
index_min(matrix, dim)

index_min(cube)
index_min(cube, dim)
```

## Examples

```cpp
[[cpp11::register]] doubles index_min1_(const int& n) {
  vec v(n, fill::randu);

  uword i = index_max(v);
  double max_val_in_v = v(i);


  mat M(n, n + 1, fill::randu);

  urowvec ii = index_max(M);
  ucolvec jj = index_max(M, 1);

  // max values in col 0 and row n
  return writable::doubles res({M(ii(0), 0), M(n, jj(n))});
}
```

# In-place dense transpose {#inplace_trans}

The `inplace_trans()` and `inplace_strans()` function return the in-place
transpose of a dense matrix. For both functions the optional `method = "lowmem"`
argument uses a low memory (and slower) algorithm for the transpose (the default
is `method = "std"`).

For real matrices:

* `inplace_trans()` returns the common transpose of the input matrix.
* `inplace_strans()` does not apply.

For complex matrices:

* `inplace_trans()` returns the Hermitian transpose (conjugate transpose) of the
  input matrix.
* `inplace_strans()` returns the transposed copy without taking the conjugate of
  the elements of the input matrix.

## Examples

```cpp
[[cpp11::register]] doubles_matrix<> inplace_trans1_(const int& n) {
  mat X(n, n, fill::randu);
  inplace_trans(X);
  return as_doubles_matrix(X);
}

[[cpp11::register]] list inplace_strans1_(const int& n) {
  cx_mat X(n, n, fill::randu);
  inplace_strans(X);
  return as_complex_matrix(X);
}
```

# Find common elements in two vectors/matrices {#intersect}

The `intersect()` function returns the common elements for two vectors or
matrices. The output is an ascending sorted vector of unique common elements.

## Examples

```cpp
[[cpp11::register]] integers intersect1_(const int& n) {
  ivec A = regspace<ivec>(n, 1);      // n, ..., 1
  ivec B = regspace<ivec>(2, n + 1);  // 2, ..., n + 1

  ivec C = intersect(A, B);  // 2, ..., n

  return as_integers(C);
}
```

# Concatenation of matrices {#join_rows}

The `join_rows()` and `join_cols()` functions concatenate matrices horizontally
and vertically, respectively. The input matrices must have the same number of
rows for `join_rows()` and the same number of columns for `join_cols()`. Both
functions accept from two to four matrices as input.

Alternatively, `join_horiz()` and `join_vert()` can be used as aliases for
`join_rows()` and `join_cols()`, respectively.

## Examples

```cpp
[[cpp11::register]] list join_rows1_(const int& n) {
  mat A(n, 1, fill::randu);
  mat B(n, 1, fill::randu);
  mat C(n, 1, fill::randu);

  mat D = join_rows(A, B, C);
  mat E = join_cols(A, B, C);

  return writable::list({A, B, C, D, E});
}
```

# Concatenation of cubes {#join_slices}

The `join_slices()` function concatenates cubes along the third dimension. For
two matrices, the input matrices must have the same number of rows and columns.
For two cubes, the input cubes must have the same number of rows and columns.
For matrix and cube, the number of rows and columns of the matrix must match the
number of rows and columns of the cube.

Usage:

```cpp
join_slices(matrix, matrix)
join_slices(cube, cube);
join_slices(matrix, cube);
join_slices(cube, matrix);
```

## Examples

```cpp
[[cpp11::register]] list join_cubes1_(const int& n) {
  cube C(n, n + 1, 3, fill::randu);
  cube D(n, n + 1, 4, fill::randu);

  cube E = join_slices(C, D);

  size_t m = C.n_slices + D.n_slices;

  writable::list res(m);

  for (size_t i = 0; i < m; ++i) {
    res[i] = as_doubles_matrix(E.slice(i));
  }

  return res;
}
```

# Kronecker tensor product {#kron}

The `kron()` function computes the Kronecker tensor product of two matrices.

## Examples

```cpp
[[cpp11::register]] doubles_matrix<> kron1_(const int& n) {
  mat A(n, n + 1, fill::randu);
  mat B(n + 1, n, fill::randu);

  mat K = kron(A, B);

  return as_doubles_matrix(K);
}
```

# Log determinant {#log_det}

The `log_det()` function computes the natural logarithm of the determinant of a
square matrix based on LU decomposition. If the matrix is not square or the
computation fails, the function throws a `std::runtime_error` exception.

Usage:

```cpp
complex val = log_det(X);
log_det(val, sign, X);
```

Form 1: `log_det(X)` returns the complex logarithm of the determinant of `X`. If
the input matrix is real, the imaginary part of the result is zero.

Form 2: `log_det(val, sign, X)` returns a bool indicating if the calculation
was successful and stores the logarithm of the determinant in the `val` and
`sign` variables such that `det(X) = sign * exp(val)`. If the computation fails,
the values of `val` and `sign` are undefined and it returns `false` without
throwing an exception.

## Examples

```cpp
[[cpp11::register]] list log_det1_(const int& n) {
  mat A(n, n, fill::randu);

  cx_double res1 = log_det(A);  // form 1

  cpp11::writable::list res2;
  res2.push_back(writable::doubles({std::real(res1)}));
  res2.push_back(writable::doubles({std::imag(res1)}));

  double val;
  double sign;
  bool ok = log_det(val, sign, A);  // form 2

  writable::list res3(3);
  res3[0] = doubles({val});
  res3[1] = doubles({sign});
  res3[2] = logicals({ok});

  writable::list res(2);
  res[0] = res2;
  res[1] = res3;

  return res;
}
```

# Log determinant of symmetric positive definite matrix {#log_det_sympd}

The `log_det_sympd()` function computes the natural logarithm of the determinant
of a symmetric positive definite matrix. If the matrix is not square or the
computation fails, a `std::runtime_error` exception is thrown.

Form 1: `log_det_sympd(X)` returns the logarithm of the determinant of `X`.

Form 2: `log_det_sympd(val, X)` returns a bool indicating if the calculation
was successful and stores the logarithm of the determinant in the `val`
variable. If the computation fails, the value of `val` is undefined and it
returns `false` without throwing an exception.

## Examples

```cpp
[[cpp11::register]] list log_det_sympd1_(const int& n) {
  mat A(n, n, fill::randu);
  A = A * A.t();  // make A symmetric positive definite

  double val = log_det_sympd(A);  // form 1

  double val2;
  bool ok = log_det_sympd(val2, A);  // form 2

  writable::list res(2);
  res[0] = doubles({val});

  writable::list res2(2);
  res2[0] = doubles({val2});
  res2[1] = logicals({ok});
  res[1] = res2;

  return res;
}
```

# Matrix logarithm {#logmat}

The `logmat()` function computes the matrix logarithm of a square matrix. If the
input matrix is not square or the computation fails, a `std::runtime_error`
exception is thrown.

Form 1: `logmat(X)` returns the matrix logarithm of `X`.

Form 2: `logmat(val, X)` returns a bool indicating if the calculation was
successful and stores the matrix logarithm in the `val` variable. If the
computation fails, the value of `val` is undefined and it returns `false` without
throwing an exception.

## Examples

```cpp
[[cpp11::register]] list logmat1_(const int& n) {
  mat A(n, n, fill::randu);
  cx_mat B = logmat(A);
  return as_complex_matrix(B);
}
```

## Caveats

* The matrix logarithm operation is generally not the same as applying the
  `log()` function to each element.
* If the input matrix is symmetric positive definite, `logmat_sympd()` is
  faster.

# Matrix logarithm of symmetric matrix {#logmat_sympd}

The `logmat_sympd()` function computes the matrix logarithm of a symmetric
positive definite matrix. If the input matrix is not square or the computation
fails, a `std::runtime_error` exception is thrown.

Form 1: `logmat_sympd(X)` returns the matrix logarithm of `X`.

Form 2: `logmat_sympd(Y, X)` returns a bool indicating if the calculation was
successful and stores the matrix logarithm in the `Y` variable. If the
computation fails, the value of `Y` is undefined and it returns `false` without
throwing an exception.

## Examples

```cpp
[[cpp11::register]] doubles_matrix<> logmat_sympd1_(const int& n) {
  mat A(n, n, fill::randu);
  mat B = A * A.t();  // make symmetric matrix
  mat C = logmat_sympd(B);
  return as_doubles_matrix(C);
}
```

# Return extremum values {#min}

The `min()` and `max()` functions return the minimum and maximum values in a
vector, matrix or cube. For a vector, the output is a scalar. For a matrix, the
output is a vector with the minimum or maximum value for each column when
`dim = 0` (default) and each row when `dim = 1`. For a cube, the output is a
cube with the minimum or maximum value for each slice's columns when `dim = 0`,
the minimum or maximum value for each slice's rows when `dim = 1`, and the
minimum or maximum value for each slice when `dim = 2`. For complex objects,
the absolute value is used to compare the elements.

Usage:

```cpp
// max() is analogous

min(vector);
min(vector1, vector2);

min(matrix);
min(matrix, dim);
min(matrix1, matrix2);

min(cube);
min(cube, dim);
min(cube1, cube2);
```

## Examples

```cpp
[[cpp11::register]] list max1_(const int& n) {
  mat M(n, n, fill::randu);

  rowvec a = max(M);
  rowvec b = max(M, 0);
  colvec c = max(M, 1);

  // element-wise maximum
  mat X(n, n, fill::randu);
  mat Y(n, n, fill::randu);
  mat Z = arma::max(X, Y);  // use arma:: prefix to distinguish from std::max()

  writable::list res(4);
  res[0] = as_doubles(a.t());
  res[1] = as_doubles(b.t());
  res[2] = as_doubles(c);
  res[3] = as_doubles_matrix(Z);

  return res;
}
```

# Return non-zero values {#nonzeros}

The `nonzeros()` function returns the non-zero values in a vector, matrix or
cube. The output is a column vector of non-zero values (`vec`). The input matrix
can be dense or sparse.

## Examples

```cpp
[[cpp11::register]] doubles nonzeros1_(const int& n) {
  mat A(n, n, fill::randu);
  A.elem(find(A < 0.5)).zeros();  // set elements less than 0.5 to zero
  vec B = nonzeros(A);
  return as_doubles(B);
}
```

## Caveats

 Caveats:

* For dense matrices/vectors, to obtain the number of non-zero elements, the
  expression `accu(X != 0)` is more efficient.
* For sparse matrices, to obtain the number of non-zero elements,
  `X.n_nonzero` is more efficient.

# Various norms of vectors and matrices {#norm}

The `norm()` function computes the p-norm of a vector or matrix. The optional
argument `p` can be `p = {1,...,n}`, `p = "inf`", `p = "-inf"`, or `p = "fro"`
for the 1,2,...,n-norms, maximum norm, minimum quasi-norm, and Frobenius norm,
respectively. The default is the 2-norm for vectors and the Frobenius norm for
matrices.

## Examples

```cpp
[[cpp11::register]] doubles norm1_(const int& n) {
  vec A(n, fill::randu);
  mat B(n, n, fill::randu);

  double a1 = norm(A, 1);
  double a2 = norm(A, 2);
  double a3 = norm(A, "inf");
  double a4 = norm(A, "-inf");
  double a5 = norm(A, "fro");

  double b1 = norm(B, 1);
  double b2 = norm(B, 2);
  double b3 = norm(B, "inf");
  double b4 = norm(B, "-inf");
  double b5 = norm(B, "fro");

  writable::doubles res({a1, a2, a3, a4, a5, b1, b2, b3, b4, b5});
  attr(res, "names") = strings({"a1", "a2", "a3", "a4", "a5",
    "b1", "b2", "b3", "b4", "b5"});
}
```

## Caveats

* The matrix 2-norm (spectral norm) is based on SVD, which is computationally
  intensive. A faster alternative is `norm2est()`.
* To obtain the vector norm of each row or column of a matrix, use `vecnorm()`.
* To obtain the zero/Hamming pseudo-norm (number of non-zero elements), use the
  expression  `accu(X != 0)`.

# Fast estimate of the matrix 2-norm {#norm2est}

The `norm2est()` function computes a fast estimate of the 2-norm of a matrix.
The function iterates until `|est1 - est2| / max(est1, est2) < tol` or the
number of iterations is equal to `max_iter`. The default values are `tol = 1e-5`
and `max_iter = 100`.

## Examples

```cpp
[[cpp11::register]] doubles norm2est1_(const int& n) {
  mat A(n, n, fill::randu);
  return doubles({norm2est(A)});
}
```

# Normalise vectors to unit p-norm {#normalise}

The `normalise()` function normalises vectors or matrices to a p-norm. The
default is the 2-norm for vectors and matrices (`p = 2`). For matrices, the
optional `dim` argument specifies the dimension along which to normalise the
matrix, with `dim = 0` normalising along columns and `dim = 1` normalising along
rows.

## Examples

```cpp
[[cpp11::register]] list normalise1_(const int& n) {
  mat A(n, n, fill::randu);

  mat B = normalise(A, 1, 0);
  mat C = normalise(A, 1, 1);

  writable::list res(2);
  res[0] = as_doubles_matrix(B);
  res[1] = as_doubles_matrix(C);

  res.attr("names") = strings({"B_norm1_cols", "C_norm1_rows"});

  return res;
}
```

# Element-wise power {#pow}

The `pow()` function computes the element-wise power of a matrix or vector. The
power argument can be a scalar, vector, or matrix.

## Examples

```cpp
[[cpp11::register]] list pow1_(const int& n) {
  mat A(n, n, fill::randu);
  mat B(n, n, fill::randu);

  mat C = pow(A, 2);
  mat D = pow(A, B);

  writable::list res(2);
  res[0] = as_doubles_matrix(C);
  res[1] = as_doubles_matrix(D);

  return res;
}
```

## Caveats

* To raise all elements to the power 2, use `square()` instead.
* For the matrix power operation, which takes into account matrix structure, use
  `powmat()`.

# Matrix power {#powmat}

The `powmat()` function computes the matrix power of a square matrix. The power
argument must be a scalar (e.g., `double` or `int`). If the input matrix is not
square, the function throws a `std::runtime_error` exception.

Usage:

```
Y = powmat(X, 2); // store a matrix
powmat(Y, X, 2); // store the matrix in Y and return true if successful
```

If the calculation fails:

* `Y = powmat(X)` throws a `std::runtime_error` exception.
* `powmat(Y, X, 2)` returns a bool set to false (exception is not thrown).

## Examples

```cpp
[[cpp11::register]] list powmat1_(const int& n) {
  mat A(n, n, fill::randu);

  mat B = powmat(A, 2);  // form 1

  mat C;
  bool ok = powmat(C, A, 2);  // form 2

  writable::list res(2);
  res[0] = as_doubles_matrix(B);

  writable::list res2(2);
  res2[0] = as_doubles_matrix(C);
  res2[1] = logicals({ok});

  res[1] = res2;

  res.attr("names") = strings({"powmat_form1", "powmat_form2"});
  res2.attr("names") = strings({"result", "status"});

  return res;
}
```

# Product of elements {#prod}

The `prod()` function computes the product of the elements in a vector or
matrix. The optional `dim` argument specifies the dimension along which to
compute the matrix product, with `dim = 0` computing the product along columns
and `dim = 1` computing the product along rows.

## Examples

```cpp
[[cpp11::register]] list prod1_(const int& n) {
  mat A(n, n, fill::randu);

  rowvec b = prod(A, 0);
  vec c = prod(A, 1);

  writable::list res(2);
  res[0] = as_doubles(b.t());
  res[1] = as_doubles(c);

  return res;
}
```

# Rank of matrix {#rank}

The `rank()` function computes the rank of a matrix based on singular values.
The optional `tolerance` argument specifies the tolerance for the singular
values. The default is `tolerance = max_rc * max_sv * epsilon`, where:

* `max_rc = max(X.n_rows, X.n_cols)`
* `max_sv = max(singular values of X)`
* `epsilon = 1 - min(singular values of X > 1)`

Usage:

```cpp
val = rank(X, tolerance); // form 1
rank(val, X, tolerance);   // form 2
```

## Examples

```cpp
[[cpp11::register]] list rank1_(const int& n) {
  mat A(n, n, fill::randu);

  int r1 = rank(A);

  uword r2;
  bool ok = rank(r2, A);

  writable::list res(2);
  res[0] = integers({r1});

  writable::list res2(2);
  res2[0] = integers({static_cast<int>(r2)});
  res2[1] = logicals({ok});

  res[1] = res2;

  res.attr("names") = strings({"rank1", "rank2"});
  res2.attr("names") = strings({"result", "status"});

  return res;
}
```

# Reciprocal condition number {#rcond}

The `rcond()` function computes the 1-norm estimate of the reciprocal condition
number of a square matrix. Values close to one indicate a well-conditioned
matrix, while values close to zero indicate a poorly conditioned matrix. If the
input matrix is not square, the function throws a `std::runtime_error`
exception.

## Examples

```cpp
[[cpp11::register]] doubles rcond1_(const int& n) {
  mat A(n, n, fill::randu);
  return doubles({rcond(A)});
}
```

## Caveat

To efficiently calculate the reciprocal condition and the matrix inverse at
the same time, use `inv()`.

# Replicate elements {#repelem}

The `repelem()` function replicates the elements of a matrix.

Usage:

```cpp
repelem(A, num_copies_per_row, num_copies_per_col)
```

The generated matrix has the following size:

* `n_rows = num_copies_per_row * A.n_rows`
* `n_cols	= num_copies_per_col * A.n_cols`

## Examples

```cpp
[[cpp11::register]] list repelem1_(const int& n) {
  mat A(n, n, fill::randu);
  mat B = repelem(A, 2, 3);

  writable::list res(2);
  res[0] = as_doubles_matrix(A);
  res[1] = as_doubles_matrix(B);

  return res;
}
```

# Replicate matrix in block-like fashion {#repmat}

The `repmat()` function replicates a matrix in a block-like fashion.

Usage:

```cpp
repmat(A, num_reps_row, num_reps_col)
```

The generated matrix has the following size:

* `n_rows = num_reps_row * A.n_rows`
* `n_cols = num_reps_col * A.n_cols`

## Examples

```cpp
[[cpp11::register]] list repmat1_(const int& n) {
  mat A(n, n, fill::randu);
  mat B = repmat(A, 2, 3);

  writable::list res(2);
  res[0] = as_doubles_matrix(A);
  res[1] = as_doubles_matrix(B);

  return res;
}
```

## Caveat

To apply a vector operation on each row or column of a matrix, it is generally
more efficient to use `.each_row()` or `.each_col()`.

# Change size while keeping elements {#reshape}

The `reshape()` function changes the size of a vector, matrix or cube while
keeping the elements in the same order.

Usage:

```cpp
reshape(vector, n_rows, n_cols)
reshape(matrix, n_rows, n_cols)

reshape(vector, size(matrix))
reshape(matrix, size(matrix))

reshape(cube, n_rows, n_cols, n_slices)
reshape(cube, size(cube))
```

## Examples

```cpp
[[cpp11::register]] list reshape1_(const int& n) {
  mat A(n, n + 1, fill::randu);
  
  mat B = reshape(A, n + 1, n);
  
  mat C(n + 4, n - 1);
  C = reshape(A, size(C));

  writable::list res(2);
  res[0] = as_doubles_matrix(B);
  res[1] = as_doubles_matrix(C);

  return res;
}
```

# Change size while keeping elements and preserving layout {#resize}

The `resize()` function changes the size of a vector, matrix or cube while
preserving the data. If the new size is larger, the new elements are set to
zero.

Usage:

```cpp
resize(vector, n_rows, n_cols)
resize(matrix, n_rows, n_cols)

resize(vector, size(matrix))
resize(matrix, size(matrix))

resize(cube, n_rows, n_cols, n_slices)
resize(cube, size(cube))
```

## Examples

```cpp
[[cpp11::register]] list resize2_(const int& n) {
  mat A(n, n + 1, fill::randu);

  mat B = resize(A, n + 1, n);

  mat C(n + 4, n - 1);
  C = resize(A, size(C));

  writable::list res(3);
  res[0] = as_doubles_matrix(A);
  res[1] = as_doubles_matrix(B);
  res[2] = as_doubles_matrix(C);

  return res;
}
```

# Reverse order of elements {#reverse}

The `reverse()` function reverses the order of elements in a vector or matrix.
The optional `dim` argument specifies the dimension along which to reverse the
matrix, with `dim = 0` reversing along columns and `dim = 1` reversing along
rows (`dim = 0` by default).

## Examples

```cpp
[[cpp11::register]] list reverse1_(const int& n) {
  mat A(n, n, fill::randu);

  mat B = reverse(A, 0);
  mat C = reverse(A, 1);

  writable::list res(3);
  res[0] = as_doubles_matrix(A);
  res[1] = as_doubles_matrix(B);
  res[2] = as_doubles_matrix(C);

  return res;
}
```

# Roots of polynomial {#roots}

The `roots()` function computes the roots of a polynomial with real or complex
coefficients. The input is a vector of coefficients, with the first element
corresponding to the highest degree term. If the computation fails, the function
throws a `std::runtime_error` exception.

Usage:

```cpp
Y = roots(X) // store the roots in Y
roots(Y, X)  // store the roots in Y and return true if successful
```

## Examples

```cpp
[[cpp11::register]] list roots1_(const int& n) {
  // y = p_1*x^n + p_2*x^(n-1) + ... + p_(n-1)*x + p_n
  // p_1, ..., p_n are random numbers
  vec y(n, 1, fill::randu);

  // note that mat and cx_mat operate directly
  // but vec and cx_vec require conv_to<...>::from()
  cx_vec z = roots(conv_to<cx_vec>::from(y));

  list res = as_complex_doubles(z);
  return res;
}
```

# Shift elements {#shift}

The `shift()` function generates a copy of a vector `V` or a matrix `M` with the
elements shifted by `N` positions in a circular manner. The `N` argument can be
positive or negative. For a matrix, the optional `dim` argument specifies the
dimension along which to shift the matrix, with `dim = 0` shifting along columns
(default) and `dim = 1` shifting along rows.

Usage:

```cpp
shift(V, N)
shift(M, N)
shift(M, N, dim)
```

## Examples

```cpp
[[cpp11::register]] list shift1_(const int& n) {
  mat A(n, n, fill::randu);
  mat B = shift(A, -1);
  mat C = shift(A, +1);

  writable::list res(3);
  res[0] = as_doubles_matrix(A);
  res[1] = as_doubles_matrix(B);
  res[2] = as_doubles_matrix(C);

  return res;
}
```

# Randomly shuffle elements {#shuffle}

<!-- shuffle( V )
shuffle( X )
shuffle( X, dim )

    For vector V, generate a copy of the vector with the elements shuffled

    For matrix X, generate a copy of the matrix with the elements shuffled in each column (dim = 0), or each row (dim = 1)

    The dim argument is optional; by default dim = 0 is used

    Examples:

        mat A(4, 5, fill::randu);
        mat B = shuffle(A);
 -->

The `shuffle()` function generates a copy of a vector `V` or matrix `M` with the
elements shuffled. For a matrix, the optional `dim` argument specifies the
dimension along which to shuffle the matrix, with `dim = 0` shuffling along
columns (default) and `dim = 1` shuffling along rows.

Usage:

```cpp
shuffle(V)
shuffle(M)
shuffle(M, dim)
```

## Examples

```cpp
[[cpp11::register]] list shuffle1_(const int& n) {
  mat A(n, n, fill::randu);
  mat B = shuffle(A);

  writable::list res(2);
  res[0] = as_doubles_matrix(A);
  res[1] = as_doubles_matrix(B);

  return res;
}
```

# Obtain dimensions of given object {#size}

The `size()` function obtains the dimensions of a matrix or cube `X`. It can
also be used to explicitly specify the dimensions of a matrix or cube.

Usage:

```cpp
size(X)
size(n_rows, n_cols)
size(n_rows, n_cols, n_slices)
```

## Examples

```cpp
[[cpp11::register]] list size1_(const int& n) {
  mat A(n, n, fill::randu);

  mat B(size(A), fill::zeros);

  mat C;
  C.randu(size(A));
  mat D = ones<mat>(size(A));

  mat E(2 * n, 2 * n, fill::ones);
  E(1, 2, size(C)) = C;  // access submatrix of E

  mat F(size(A) + size(E), fill::randu);

  mat G(size(A) * 2, fill::randu);

  writable::list res(7);

  res[0] = as_doubles_matrix(A);
  res[1] = as_doubles_matrix(B);
  res[2] = as_doubles_matrix(C);
  res[3] = as_doubles_matrix(D);
  res[4] = as_doubles_matrix(E);
  res[5] = as_doubles_matrix(F);
  res[6] = as_doubles_matrix(G);

  return res;
}
```

# Sort elements {#sort}

The `sort()` function returns a sorted version of a vector `V` or matrix `M`.
For a matrix, the optional `dim` argument specifies the dimension along which to
sort the matrix, with `dim = 0` sorting along columns (default) and `dim = 1`
sorting along rows. The optional `sort_direction` argument specifies the sorting
direction, with `sort_direction = "ascend"` (default) sorting in ascending order
and `sort_direction = "descend"` sorting in descending order.

Usage:

```cpp
sort(V)
sort(V, sort_direction)

sort(M)
sort(M, sort_direction)
sort(M, sort_direction, dim)
```

## Examples

```cpp
[[cpp11::register]] list sort1_(const int& n) {
  mat A(n, n, fill::randu);
  mat B = sort(A);
  mat C = sort(A, "descend");
  mat D = sort(A, "ascend", 1);
  mat E = sort(A, "descend", 1);

  writable::list res(5);
  res[0] = as_doubles_matrix(A);
  res[1] = as_doubles_matrix(B);
  res[2] = as_doubles_matrix(C);
  res[3] = as_doubles_matrix(D);
  res[4] = as_doubles_matrix(E);

  return res;
}
```

# Vector describing sorted order of elements {#sort_index}

The `sort_index()` function returns a vector describing the sorted order of the
elements of a vector `V` or matrix `M`. The optional `sort_direction` argument
specifies the sorting direction, with `sort_direction = "ascend"` (default)
sorting in ascending order and `sort_direction = "descend"` sorting in
descending order.

Usage:

```cpp
sort_index(V)
sort_index(V, sort_direction)

sort_index(M)
sort_index(M, sort_direction)
```

## Examples

```cpp
[[cpp11::register]] list sort_index1_(const int& n) {
  mat A(n, n, fill::randu);
  uvec B = sort_index(A);
  uvec C = sort_index(A, "descend");

  writable::list res(3);
  res[0] = as_doubles_matrix(A);
  res[1] = as_integers(B);
  res[2] = as_integers(C);

  return res;
}
```

# Generate a sparse matrix with diagonals specified by column vectors {#spdiags}

The `spdiags()` function generates a sparse matrix with diagonals specified by
column vectors from an input matrix and a vector to indicate the diagonals.

Usage:

```cpp
spdiags(matrix, vector, number_of_rows, number_of_columns);
```

Each element in the input vector specifies diagonal `k`, where:

* `k = 0` is the main diagonal
* `k > 0` is above the main diagonal
* `k < 0` is below the main diagonal

## Examples

```cpp
[[cpp11::register]] doubles_matrix<> spdiags1_(const int& n) {
  mat V(n, n, fill::randu);
  ivec D = {0, -1};
  sp_mat X = spdiags(V, D, n, n); // lower triangular matrix
  return as_doubles_matrix(X);
}
```

# Square root of matrix {#sqrtmat}

The `sqrtmat()` function computes the complex square root of a general square
matrix. If the input matrix is not square, the function throws an error. If the
matrix appears to be singular, an approximate square root is attempted.

Usage:

```cpp
B = sqrtmat(A)
sqrtmat(B, A)
```

## Examples

```cpp
[[cpp11::register]] list sqrtmat1_(const int& n) {
  mat A(n, n, fill::randu);
  
  cx_mat B = sqrtmat(A);

  cx_mat C;
  bool ok = sqrtmat(C, A);

  writable::list res(4);

  res[0] = as_doubles_matrix(A);
  res[1] = as_complex_matrix(B);
  res[2] = as_complex_matrix(C);
  res[3] = logicals({ok});

  return res;
}
```

# Square root of symmetric matrix {#sqrtmat_sympd}

The `sqrtmat_sympd()` function computes the square root of a symmetric positive
definite matrix. If the input matrix is not square or the computation fails, the
function throws an error.

Usage:

```cpp
B = sqrtmat_sympd(A)
sqrtmat_sympd(B, A)
```

## Examples

```cpp
[[cpp11::register]] doubles_matrix<> sqrtmat_sympd1_(const int& n) {
  mat A(n, n, fill::randu);
  A = A * A.t();  // make A symmetric positive definite

  mat B = sqrtmat_sympd(A);

  return as_doubles_matrix(B);
}
```

# Sum of elements {#sum}

The `sum()` function computes the sum of the elements in a vector, matrix or
cube. For a matrix, the optional `dim` argument specifies the dimension along
which to compute the sum, with `dim = 0` computing the sum along columns and
`dim = 1` computing the sum along rows. For a cube, the optional `dim` argument
specifies the dimension along which to compute the sum, with `dim = 0` computing
the sum along columns, `dim = 1` computing the sum along rows, and `dim = 2`
computing the sum along slices.

Usage:

```cpp
sum(vector)

sum(matrix)
sum(matrix, dim)

sum(cube)
sum(cube, dim)
```

## Examples

```cpp
[[cpp11::register]] list sum2_(const int& n) {
  mat A(n, n, fill::randu);

  vec a = sum(A, 1);
  vec b = sum(A, 0).t();
  double c = accu(A);  // overall sum

  writable::list res(3);
  res[0] = as_doubles(a);
  res[1] = as_doubles(b);
  res[2] = doubles({c});

  return res;
}
```

# Convert subscripts to linear index {#sub2ind}

The `sub2ind()` function converts subscripts to a linear index. If a subscript
is out of range, the function returns an error.

Usage:

```cpp
sub2ind(size(matrix), row, col)
sub2ind(size(matrix), matrix_of_subscripts)

sub2ind(size(cube), row, col, slice)
sub2ind(size(cube), matrix_of_subscripts)
```

## Examples

```cpp
[[cpp11::register]] integers sub2ind1_(const int& n) {
  mat M(n, n, fill::randu);

  uword i = sub2ind(size(M), n - 1, n - 1);

  return integers({static_cast<int>(i)});
}
```

# Generate symmetric matrix from given matrix {#symmatu-symmatl}

The `symmatu()` function generates a symmetric matrix from a square matrix `A`
by reflecting the upper triangle to the lower triangle. The `symmatl()` function
generates a symmetric matrix from a square matrix `A` by reflecting the lower
triangle to the upper triangle. If `A` is a complex matrix, the reflection uses
the complex conjugate of the elements. To disable the complex conjugate, set
`do_conj` to `false`. If `A` is non-square, an error is thrown.

Usage:

```cpp
symmatu(A)
symmatu(A, do_conj)

symmatl(A)
symmatl(A, do_conj)
```

## Examples

```cpp
[[cpp11::register]] doubles_matrix<> symmatu1_(const int& n) {
  mat A(n, n, fill::randu);
  mat B = symmatu(A);
  return as_doubles_matrix(B);
}
```

# Sum of diagonal elements {#trace}

The `trace()` function computes the sum of the elements on the main diagonal of
a matrix. If the input matrix is not square, an error is thrown.

Usage:

```cpp
trace(X)
```

## Examples

```cpp
[[cpp11::register]] doubles trace1_(const int& n) {
  mat A(n, n, fill::randu);
  return doubles({trace(A)});
}
```

# Transpose of matrix {#trans-strans}

The `trans()` function transposes a matrix. For a real matrix,
`trans()` provides a transposed copy of the matrix. For a complex matrix,
`trans()` provides a Hermitian (conjugate) transposed copy, where the signs of
the imaginary components are flipped. The `strans()` function provides a simple
transposed copy, where the signs of the imaginary components are not flipped.

Usage:

```
trans(A)
strans(A)
```

## Examples

```cpp
[[cpp11::register]] list trans1_(const int& n) {
  mat A(n, n, fill::randu);
  
  mat B = trans(A);
  mat C = A.t();  // same as trans(A)

  writable::list res(2);

  res[0] = as_doubles_matrix(A);
  res[1] = as_doubles_matrix(C);

  return res;
}
```

# Trapezoidal numerical integration {#trapz}

The `trapz()` function computes the trapezoidal integral of a vector `Y` with
respect to spacing in a vector `X`. The optional `dim` argument specifies the
dimension along which to compute the trapezoidal integral, with `dim = 0`
computing the integral along columns and `dim = 1` computing the integral along
rows.

Usage:

```cpp
trapz(X, Y)
trapz(X, Y, dim)

trapz(Y)
trapz(Y, dim)
```

## Examples

```cpp
[[cpp11::register]] doubles_matrix<> trapz1_(n) {
  vec X = linspace<vec>(0, datum::pi, n);
  vec Y = sin(X);
  
  mat Z = trapz(X,Y);

  return as_doubles_matrix(Z);
}
```

# Copy upper/lower triangular part {#trimatu-trimatl}

The `trimatu()` function creates a new matrix by copying the upper triangular
part from a square matrix `A` and setting the remaining elements to zero. The
`trimatl()` function creates a new matrix by copying the lower triangular part
from a square matrix `A` and setting the remaining elements to zero. The
optional `k` argument specifies the diagonal (`k = 0` by default, which sets
the main diagonal). For `k > 0`, the `k`-th upper-diagonal is used (above the
main diagonal, towards the top-right corner). For `k < 0`, the `k`-th
lower-diagonal is used (below the main diagonal, towards the bottom-left
corner).

Usage:

```cpp
trimatu(A)
trimatu(A, k)

trimatl(A)
trimatl(A, k)
```

## Examples

```cpp
[[cpp11::register]] doubles_matrix<> trimatu1_(const int& n) {
  mat A(n, n, fill::randu);
  mat B = trimatu(A);
  return as_doubles_matrix(B);
}
```

# Obtain indices of upper/lower triangular part {#trimatu_ind-trimatl_ind}

The `trimatu_ind()` function returns a column vector containing the indices of
elements that form the upper triangular part of a matrix `A`. The
`trimatl_ind()` function returns a column vector containing the indices of
elements that form the lower triangular part of a matrix `A`. The optional `k`
argument specifies the diagonal (`k = 0` by default, which sets the main
diagonal). For `k > 0`, the `k`-th upper-diagonal is used (above the main
diagonal, towards the top-right corner). For `k < 0`, the `k`-th lower-diagonal
is used (below the main diagonal, towards the bottom-left corner).

Usage:

```cpp
trimatu_ind(size(A))
trimatu_ind(size(A), k)

trimatl_ind(size(A))
trimatl_ind(size(A), k)
```

## Examples

```cpp
[[cpp11::register]] integers trimatu_ind1_(const int& n) {
  mat A(n, n, fill::randu);
  uvec B = trimatu_ind(size(A));
  return as_integers(B);
}
```

# Return unique elements {#unique}

The `unique()` function returns the unique elements of a vector or matrix `A`,
sorted in ascending order. If `A` is a vector, the output is also a vector with
the same orientation (row or column) as `A`. If `A` is a matrix, the output is
always a column vector.

Usage:

```cpp
unique(A)
```

## Examples

```cpp
[[cpp11::register]] doubles unique1_(const int& n) {
  mat A(n, n, fill::randu);
  A(0, 0) = A(1, 1)
  vec B = unique(A);
  return as_doubles(B);
}
```

# Obtain vector norm of each row or column of a matrix {#vecnorm}

The `vecnorm()` function computes the p-norm of each column vector (when
`dim = 0`) or row vector (when `dim = 1`) of a matrix `X`. The optional `p`
argument specifies the norm to compute, with `p = 2` (default) computing the
2-norm, `p = 1` computing the 1-norm, `p = "inf"` computing the maximum norm,
and `p = "-inf"` computing the minimum quasi-norm.

Usage:

```cpp
vecnorm(X)
vecnorm(X, p)
vecnorm(X, p, dim)
```

## Examples

```cpp
[[cpp11::register]] list vecnorm1_(const int& n) {
  mat A(n, n, fill::randu);

  colvec a = vecnorm(A, 2).t();
  colvec b = vecnorm(A, "inf", 1);

  writable::list res(2);
  res[0] = as_doubles(a);
  res[1] = as_doubles(b);

  return res;
}
```

# Flatten matrix into vector {#vectorise}

The `vectorise()` function generates a flattened version of a matrix `M` or cube
`Q`. The optional `dim` argument specifies the dimension along which to flatten
the matrix, with `dim = 0` flattening column-wise (default) and `dim = 1`
flattening row-wise.

Usage:

```cpp
vectorise(M)
vectorise(M, dim)

vectorise(Q)
```

## Examples

```cpp
[[cpp11::register]] doubles vectorise1_(const int& n) {
  mat A(n, n, fill::randu);
  vec B = vectorise(A);
  return as_doubles(B);
}
```

# Miscellaneous element-wise functions: exp, log, sqrt, round, sign, and others {#misc-functions}

Miscellaneous element-wise functions include:

| Function            | Description               |
|---------------------|---------------------------|
| `exp()`             | Base-e exponential: `e^x` |
| `exp2()`            | Base-2 exponential: `2^x` |
| `exp10()`           | Base-10 exponential: `10^x` |
| `expm1()`           | Compute `exp(A)-1` accurately for values of `A` close to zero (only for float and double elements) |
| `trunc_exp()`       | Base-e exponential, truncated to avoid infinity (only for float and double elements) |
| `log()`             | Natural log: `loge(x)`  |
| `log2()`            | Base-2 log: `log2(x)`   |
| `log10()`           | Base-10 log: `log10(x)` |
| `log1p()`           | Compute `log(1+A)` accurately for values of `A` close to zero (only for float and double elements) |
| `trunc_log()`       | Natural log, truncated to avoid +/-infinity (only for float and double elements) |
| `square()`          | Square: `x^2`           |
| `sqrt()`            | Square root: `x^(1.2)`  |
| `cbrt()`            | Cube root: `x^(1/3)`    |
| `floor()`           | Largest integral value that is not greater than the input value |
| `ceil()`            | Smallest integral value that is not less than the input value |
| `round()`           | Round to nearest integer, with halfway cases rounded away from zero |
| `trunc()`           | Round to nearest integer, towards zero |
| `erf()`             | Error function (only for float and double elements) |
| `erfc()`            | Complementary error function (only for float and double elements) |
| `tgamma()`          | Gamma function (only for float and double elements) |
| `lgamma()`          | Natural log of the absolute value of gamma function (only for float and double elements) |
| `sign()`            | Signum function; for each element `a` in `A`, the corresponding element `b` in `B` is: `-1` if `a < 0`, `0` if `a = 0`, `+1` if `a > 0`. If `a` is complex and non-zero, then `b = a / abs(a)` |

## Caveats

All of the above functions are applied element-wise, where each element is
treated independently. `expmat()`, `logmat()`, `sqrtmat()`, and `powmat()` take
into account matrix structure.

## Examples

```cpp
[[cpp11::register]] list misc1_(const int& n) {
  mat A(n, n, fill::randu);
  mat B = exp(A);
  mat C = log(A);
  mat D = sqrt(A);
  mat E = round(A);
  mat F = sign(A);

  writable::list res(6);
  res[0] = as_doubles_matrix(A);
  res[1] = as_doubles_matrix(B);
  res[2] = as_doubles_matrix(C);
  res[3] = as_doubles_matrix(D);
  res[4] = as_doubles_matrix(E);
  res[5] = as_doubles_matrix(F);

  return res;
}
```

# Trigonometric element-wise functions: cos, sin, tan, and others {#trig-functions}


<!-- trigonometric element-wise functions (cos, sin, tan, ...)

    For single argument functions, B = trig_fn(A), where trig_fn is applied to each element in A, with trig_fn as one of:
        cos, acos, cosh, acosh
        sin, asin, sinh, asinh
        tan, atan, tanh, atanh
        sinc, defined as sinc(x) = sin(πx) / (πx) for x ≠ 0, and sinc(x) = 1 for x = 0

    For dual argument functions, apply the function to each tuple of two corresponding elements in X and Y:
        Z = atan2(Y, X)
        Z = hypot(X, Y)

    Examples:

        mat X(5, 5, fill::randu);
        mat Y = cos(X);
 -->

Trigonometric element-wise functions include:

| Function            | Description               |
|---------------------|---------------------------|
| `cos()`             | Cosine: `cos(x)`         |
| `acos()`            | Inverse cosine: `arccos(x)` |
| `cosh()`            | Hyperbolic cosine: `cosh(x)` |
| `acosh()`           | Inverse hyperbolic cosine: `arccosh(x)` |
| `sin()`             | Sine: `sin(x)`           |
| `asin()`            | Inverse sine: `arcsin(x)` |
| `sinh()`            | Hyperbolic sine: `sinh(x)` |
| `asinh()`           | Inverse hyperbolic sine: `arcsinh(x)` |
| `tan()`             | Tangent: `tan(x)`         |
| `atan()`            | Inverse tangent: `arctan(x)` |
| `tanh()`            | Hyperbolic tangent: `tanh(x)` |
| `atanh()`           | Inverse hyperbolic tangent: `arctanh(x)` |
| `sinc()`            | Sinc function: `sinc(x) = sin(datum::pi * x) / (datum::pi * x)` for `x != 0`, and `sinc(x) = 1` for `x = 0` |
| `atan2()`           | Two-argument arctangent: `atan2(y, x)` |
| `hypot()`           | Hypotenuse: `hypot(x, y)` |

## Caveats

All of the above functions are applied element-wise, where each element is
treated independently.

## Examples

```cpp
[[cpp11::register]] list trig1_(const int& n) {
  mat A(n, n, fill::randu);
  mat B = cos(A);
  mat C = sin(A);
  mat D = tan(A);
  mat E = atan2(C, B);
  mat F = hypot(B, C);

  writable::list res(6);
  res[0] = as_doubles_matrix(A);
  res[1] = as_doubles_matrix(B);
  res[2] = as_doubles_matrix(C);
  res[3] = as_doubles_matrix(D);
  res[4] = as_doubles_matrix(E);
  res[5] = as_doubles_matrix(F);

  return res;
}
```