Absolute value

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

Usage:

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.

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

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

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

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 Ahas 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  ⎦

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

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:

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

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

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:

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

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

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:

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)

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

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:

real_object_type Y = arg(X);

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

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

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

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.

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

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.

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

Obtain complex conjugate of each element

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

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

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:

conv_to<type>::from(X) 

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

Cross product

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

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

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:

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

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

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:

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

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

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)

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

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:

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

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

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

Usage:

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

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

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:

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

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

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:

diff(vector);
diff(vector, k);

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

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

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.

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

Obtain distance of each element to next largest floating point representation

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.

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

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.

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

Matrix exponential of symmetric matrix

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.

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

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:

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.

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

Find indices of finite elements

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

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

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

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

Find indices of NaN elements

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

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

Find indices of unique elements

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

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

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.

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

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

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

Convert linear index to subscripts

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:

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)

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

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:

// index_max is analogous

index_min(vector)

index_min(matrix)
index_min(matrix, dim)

index_min(cube)
index_min(cube, dim)

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

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.

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

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

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

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.

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

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:

join_slices(matrix, matrix)
join_slices(cube, cube);
join_slices(matrix, cube);
join_slices(cube, matrix);

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

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

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

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:

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.

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

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.

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

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.

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

Matrix logarithm of symmetric matrix

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.

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

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:

// max() is analogous

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

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

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

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

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.

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

Various norms of vectors and matrices

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.

[[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"});
}

Fast estimate of the matrix 2-norm

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.

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

Normalise vectors to unit p-norm

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.

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

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

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

Matrix power

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

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

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.

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

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:

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

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

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.

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

Replicate elements

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

Usage:

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

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

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

Usage:

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

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

Change size while keeping elements

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

Usage:

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

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

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:

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

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

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

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

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:

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

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

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:

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

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

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:

shuffle(V)
shuffle(M)
shuffle(M, dim)

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

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:

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

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

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:

sort(V)
sort(V, sort_direction)

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

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

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:

sort_index(V)
sort_index(V, sort_direction)

sort_index(M)
sort_index(M, sort_direction)

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

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:

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

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

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:

B = sqrtmat(A)
sqrtmat(B, A)

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

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:

B = sqrtmat_sympd(A)
sqrtmat_sympd(B, A)

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

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:

sum(vector)

sum(matrix)
sum(matrix, dim)

sum(cube)
sum(cube, dim)

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

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

Usage:

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

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

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

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:

symmatu(A)
symmatu(A, do_conj)

symmatl(A)
symmatl(A, do_conj)

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

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:

trace(X)

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

Transpose of matrix

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)

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

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:

trapz(X, Y)
trapz(X, Y, dim)

trapz(Y)
trapz(Y, dim)

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

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:

trimatu(A)
trimatu(A, k)

trimatl(A)
trimatl(A, k)

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

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:

trimatu_ind(size(A))
trimatu_ind(size(A), k)

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

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

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:

unique(A)

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

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:

vecnorm(X)
vecnorm(X, p)
vecnorm(X, p, dim)

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

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:

vectorise(M)
vectorise(M, dim)

vectorise(Q)

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

Miscellaneous element-wise functions include:

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

Trigonometric element-wise functions include:

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