Create a matrix using data from writable auxiliary memory. By default the matrix allocates its own memory and copies data from the auxiliary memory (for safety). However, if copy_aux_mem is set to false, the matrix will instead directly use the auxiliary memory (ie. no copying). This is faster, but can be dangerous unless you know what you're doing!

The strict variable comes into effect only if copy_aux_mem is set to false (ie. the matrix is directly using auxiliary memory). If strict is set to true, the matrix will be bound to the auxiliary memory for its lifetime; the number of elements in the matrix can't be changed (directly or indirectly). If strict is set to false, the matrix will not be bound to the auxiliary memory for its lifetime, ie., the size of the matrix can be changed. If the requested number of elements is different to the size of the auxiliary memory, new memory will be allocated and the auxiliary memory will no longer be used.

Create a matrix by copying data from read-only auxiliary memory.

Create a fixed size matrix, with the size specified via template arguments. Memory for the matrix is allocated at compile time. This is generally faster than dynamic memory allocation, but the size of the matrix can't be changed afterwards (directly or indirectly).

For convenience, there are several pre-defined typedefs for each matrix type (where the types are: umat, imat, fmat, mat, cx_fmat, cx_mat). The typedefs specify a square matrix size, ranging from 2x2 to 9x9. The typedefs were defined by simply appending a two digit form of the size to the matrix type -- for example, mat33 is equivalent to mat::fixed<3,3>, while cx_mat44 is equivalent to cx_mat::fixed<4,4>.

Create a fixed size matrix, with the size specified via template arguments, and copying data from auxiliary memory.

mat A = randu<mat>(5,5);
double x = A(1,2);

mat B = A + A;
mat C = A * B;
mat D = A % B;

cx_mat X(A,B);

B.zeros();
B.set_size(10,10);
B.ones(5,6);

mat E(4,5, fill::zeros);

//
// fixed size matrices:

mat::fixed<5,6> F;
F.ones();

mat44 G;
G.randn();

cout << mat22().randu() << endl;

//
// constructing matrices from
// auxiliary (external) memory:

double aux_mem[24];
mat H(aux_mem, 4, 6, false);

mat A(5,5);  A = 123.0;

mat A(5,5);  A.fill(123.0);

Create a column vector using data from writable auxiliary memory. By default the vector allocates its own memory and copies data from the auxiliary memory (for safety). However, if copy_aux_mem is set to false, the vector will instead directly use the auxiliary memory (ie. no copying). This is faster, but can be dangerous unless you know what you're doing!

The strict variable comes into effect only if copy_aux_mem is set to false (ie. the vector is directly using auxiliary memory). If strict is set to true, the vector will be bound to the auxiliary memory for its lifetime; the number of elements in the vector can't be changed (directly or indirectly). If strict is set to false, the vector will not be bound to the auxiliary memory for its lifetime, ie., the vector's size can be changed. If the requested number of elements is different to the size of the auxiliary memory, new memory will be allocated and the auxiliary memory will no longer be used.

Create a column vector by copying data from read-only auxiliary memory.

Create a fixed size column vector, with the size specified via the template argument. Memory for the vector is allocated at compile time. This is generally faster than dynamic memory allocation, but the size of the vector can't be changed afterwards (directly or indirectly).

For convenience, there are several pre-defined typedefs for each vector type (where the types are: uvec, ivec, fvec, vec, cx_fvec, cx_vec as well as the corresponding colvec versions). The pre-defined typedefs specify vector sizes ranging from 2 to 9. The typedefs were defined by simply appending a single digit form of the size to the vector type -- for example, vec3 is equivalent to vec::fixed<3>, while cx_vec4 is equivalent to cx_vec::fixed<4>.

Create a fixed size column vector, with the size specified via the template argument, and copying data from auxiliary memory.

vec x(10);
vec y = zeros<vec>(10);

mat A = randu<mat>(10,10);
vec z = A.col(5); // extract a column vector

vec a(5);  a = 123.0;

vec a(5);  a.fill(123.0);

Create a row vector using data from writable auxiliary memory. By default the vector allocates its own memory and copies data from the auxiliary memory (for safety). However, if copy_aux_mem is set to false, the vector will instead directly use the auxiliary memory (ie. no copying). This is faster, but can be dangerous unless you know what you're doing!

The strict variable comes into effect only if copy_aux_mem is set to false (ie. the vector is directly using auxiliary memory). If strict is set to true, the vector will be bound to the auxiliary memory for its lifetime; the number of elements in the vector can't be changed (directly or indirectly). If strict is set to false, the vector will not be bound to the auxiliary memory for its lifetime, ie., the vector's size can be changed. If the requested number of elements is different to the size of the auxiliary memory, new memory will be allocated and the auxiliary memory will no longer be used.

Create a row vector by copying data from read-only auxiliary memory.

Create a fixed size row vector, with the size specified via the template argument. Memory for the vector is allocated at compile time. This is generally faster than dynamic memory allocation, but the size of the vector can't be changed afterwards (directly or indirectly).

For convenience, there are several pre-defined typedefs for each vector type (where the types are: urowvec, irowvec, frowvec, rowvec, cx_frowvec, cx_rowvec). The pre-defined typedefs specify vector sizes ranging from 2 to 9. The typedefs were defined by simply appending a single digit form of the size to the vector type -- for example, rowvec3 is equivalent to rowvec::fixed<3>, while cx_rowvec4 is equivalent to cx_rowvec::fixed<4>.

Create a fixed size row vector, with the size specified via the template argument, and copying data from auxiliary memory.

rowvec x(10);
rowvec y = zeros<rowvec>(10);

mat    A = randu<mat>(10,10);
rowvec z = A.row(5); // extract a row vector

rowvec r(5);  r = 123.0;

rowvec r(5);  r.fill(123.0);

Create a fixed size cube, with the size specified via template arguments. Memory for the cube is allocated at compile time. This is generally faster than dynamic memory allocation, but the size of the cube can't be changed afterwards (directly or indirectly).

Create a cube using data from writable auxiliary memory. By default the cube allocates its own memory and copies data from the auxiliary memory (for safety). However, if copy_aux_mem is set to false, the cube will instead directly use the auxiliary memory (ie. no copying). This is faster, but can be dangerous unless you know what you're doing!

The strict variable comes into effect only if copy_aux_mem is set to false (ie. the cube is directly using auxiliary memory). If strict is set to true, the cube will be bound to the auxiliary memory for its lifetime; the number of elements in the cube can't be changed (directly or indirectly). If strict is set to false, the cube will not be bound to the auxiliary memory for its lifetime, ie., the size of the cube can be changed. If the requested number of elements is different to the size of the auxiliary memory, new memory will be allocated and the auxiliary memory will no longer be used.

Create a cube by copying data from read-only auxiliary memory.

cube x(1,2,3);
cube y = randu<cube>(4,5,6);

mat A = y.slice(1);  // extract a slice from the cube
                     // (each slice is a matrix)

mat B = randu<mat>(4,5);
y.slice(2) = B;     // set a slice in the cube

cube q = y + y;     // cube addition
cube r = y % y;     // element-wise cube multiplication

cube::fixed<4,5,6> f;
f.ones();

// create a field containing vectors
field<vec> F(3,2);

// each vector in the field can have an arbitrary size
F(0,0) = vec(5);
F(1,1) = randu<vec>(6);
F(2,0).set_size(7);

// access element 1 of vector stored at 2,0
double x = F(2,0)(1);

// copy a row of vectors
F.row(0) = F.row(2);

// extract a row of vectors from F
field<vec> G = F.row(1);

// print the field to the standard output
G.print("G =");

// save the field to a binary file
G.save("vec_field");

• Using batch insertion constructors is generally much faster than consecutively inserting values using element access operators


• For form 1, rowind is a dense column vector of type umat or uvec containing the row indices of the values to be inserted, and colptr is a dense column vector of type umat or uvec containing indices of values corresponding to the start of new columns; the vectors correspond to the arrays used by the compressed sparse column format


• For forms 2, 3, 4, locations is a dense matrix of type umat, with a size of 2 x N, where N is the number of values to be inserted; the location of the i-th element is specified by the contents of the i-th column of the locations matrix, where the row is in locations(0,i), and the column is in locations(1,i)


• For all forms, values is a dense column vector containing the values to be inserted; it must have the same element type as the sparse matrix. For forms 2 and 3, the value in values[i] will be inserted at the location specified by the i-th column of the locations matrix.


• For form 4, add_values is either true or false; when set to true, identical locations are allowed, and the values at identical locations are added


• The size of the constructed matrix is either manually specified via n_rows and n_cols (forms 1, 3, 4), or automatically determined from the maximal locations in the locations matrix (form 2)


• If sort_locations is set to false, the locations matrix is assumed to contain locations that are already sorted according to column-major ordering


• If check_for_zeros is set to false, the values vector is assumed to contain no zero values

sp_mat A(5,6);
sp_mat B(6,5);

A(0,0) = 1;
A(1,0) = 2;

B(0,0) = 3;
B(0,1) = 4;

sp_mat C = 2*B;

sp_mat D = A*C;


// batch insertion of two values at (5, 6) and (9, 9)
umat locations;
locations << 5 << 9 << endr
          << 6 << 9 << endr;

vec values;
values << 1.5 << 3.2 << endr;

sp_mat X(locations, values);

mat X(4,5);
cout << "X has " << X.n_cols << " columns" << endl;

mat A = randu<mat>(5,5);

double* mem = A.colptr(2);

mat A = randu<mat>(5,6);
mat B;
B.copy_size(A);

cout << B.n_rows << endl;
cout << B.n_cols << endl;

mat X = randu<mat>(5,5);

vec a = X.diag();
vec b = X.diag(1);
vec c = X.diag(-2);

X.diag()  = randu<vec>(5);
X.diag() += 6;

mat X = ones<mat>(6,5);
vec v = linspace<vec>(10,15,6);

// add v to each column in X
X.each_col() += v;

// subtract v from columns 0 through to 3 in X
X.cols(0,3).each_col() -= v;

uvec indices(2);
indices(0) = 2;
indices(1) = 4;

// copy v to columns 2 and 4 in X
X.each_col(indices) = v;

mat A = randu<mat>(10,10);
A(9,9) = 123.0;
double x = A.at(9,9);
double y = A[99];

vec p = randu<vec>(10,1);
p(9) = 123.0;
double z = p[9];

mat A;

A << 1 << 2 << 3 << endr
  << 4 << 5 << 6 << endr;

mat B = { 1, 2, 3, 4, 5, 6 };  // C++11 only
B.reshape(2,3);

cx_mat A( randu<mat>(4,4), randu<mat>(4,4) );

real(A).eval().save("A_real.dat", raw_ascii);
imag(A).eval().save("A_imag.dat", raw_ascii);

mat A(5,5);
A.eye();

mat B;
B.eye(5,5);

mat A(5,5);
A.fill(123.0);

mat A = randu<mat>(4,4);
mat B = randu<mat>(4,1);

mat X = A.i();

mat Y = (A+A).i();

mat A = randu<mat>(4,5);

cout << A.in_range(0,0) << endl;  // true
cout << A.in_range(3,4) << endl;  // true
cout << A.in_range(4,5) << endl;  // false

mat A = randu<mat>(5,5);
cout << A.is_empty() << endl;

A.reset();
cout << A.is_empty() << endl;

mat A = randu<mat>(5,5);
mat B = randu<mat>(5,5);

B(1,1) = datum::inf;

cout << A.is_finite() << endl;
cout << B.is_finite() << endl;

mat A = randu<mat>(5,5);
mat B = randu<mat>(6,7);

cout << A.is_square() << endl;
cout << B.is_square() << endl;

mat A = randu<mat>(1,5);
mat B = randu<mat>(5,1);
mat C = randu<mat>(5,5);

cout << A.is_vec() << endl;
cout << B.is_vec() << endl;
cout << C.is_vec() << endl;

// C++11 only example
// need to include <random>

std::mt19937 engine;  // Mersenne twister random number engine

std::uniform_real_distribution<double> distr(0.0, 1.0);
  
mat A(4,5);
  
A.imbue( [&]() { return distr(engine); } );

mat A = randu<mat>(5,10);
mat B = ones<mat>(5,2);

// at column 2, insert a copy of B;
// A will now have 12 columns
A.insert_cols(2, B);

// at column 1, insert 5 zeroed columns;
// B will now have 7 columns
B.insert_cols(1, 5);

mat X = randu<mat>(5,5);


mat::iterator a = X.begin();
mat::iterator b = X.end();

for(mat::iterator i=a; i!=b; ++i)
  {
  cout << *i << endl;
  }


mat::col_iterator c = X.begin_col(1);  // start of column 1
mat::col_iterator d = X.end_col(3);    // end of column 3

for(mat::col_iterator i=c; i!=d; ++i)
  {
  cout << *i << endl;
  (*i) = 123.0;
  }

cube X = randu<cube>(2,3,4);


cube::iterator a = X.begin();
cube::iterator b = X.end();

for(cube::iterator i=a; i!=b; ++i)
  {
  cout << *i << endl;
  }


cube::slice_iterator c = X.begin_slice(1);  // start of slice 1
cube::slice_iterator d = X.end_slice(2);    // end of slice 2

for(cube::slice_iterator i=c; i!=d; ++i)
  {
  cout << *i << endl;
  (*i) = 123.0;
  }

      mat A = randu<mat>(5,5);
const mat B = randu<mat>(5,5);

      double* A_mem = A.memptr();
const double* B_mem = B.memptr();

vec v = randu<vec>(10);

cout << "min value is " << v.min() << endl;


uword  index;
double min_val = v.min(index);

cout << "index of min value is " << index << endl;


mat A = randu<mat>(5,5);

uword  row;
uword  col;
double min_val2 = A.max(row,col);

cout << "max value is at " << row << ',' << col << endl;

mat A = randu<mat>(5,10);
A.ones();      // sets all elements to one
A.ones(10,20); // sets the size to 10 rows and 20 columns
               // followed by setting all elements to one

mat A = randu<mat>(5,10);
mat B = randu<mat>(5,10);
mat C = randu<mat>(10,5);

mat P = A + B;
mat Q = A - B;
mat R = -B;
mat S = A / 123.0;
mat T = A % B;
mat U = A * C;

// V is constructed without temporaries
mat V = A + B + A + B;

imat AA = "1 2 3; 4 5 6; 7 8 9;";
imat BB = "3 2 1; 6 5 4; 9 8 7;";

// compare elements
umat ZZ = (AA >= BB);

mat A = randu<mat>(5,5);
mat B = randu<mat>(6,6);

A.print();

// print a transposed version of A
A.t().print();

// "B:" is the optional header line
B.print("B:");

cout << A << endl;

cout << "B:" << endl;
cout << B << endl;

mat A = randu<mat>(5,5);

cout.precision(11);
cout.setf(ios::fixed);

A.raw_print(cout, "A =");

mat A(4,5);
A.randu();

mat B;
B.randu(6,7);

arma_rng::set_seed_random();  // set the seed to a random value

mat A = randu<mat>(5, 5);
A.reset();

mat A = randu<mat>(4,5);

A.reshape(5,4);

mat A = randu<mat>(4,5);
A.resize(7,6);

mat A = randu<mat>(5,5);

A.save("A1.mat");  // default save format is arma_binary
A.save("A2.mat", arma_ascii);

mat B;
// automatically detect format type
B.load("A1.mat");

mat C;
// force loading in the arma_ascii format
C.load("A2.mat", arma_ascii);


// example of saving/loading using a stream
std::stringstream s;
A.save(s);

mat D;
D.load(s);


// example of testing for success
mat E;
bool status = E.load("A2.mat");

if(status == true)
  {
  cout << "loaded okay" << endl;
  }
else
  {
  cout << "problem with loading" << endl;
  }

   mat A = randu<mat>(4,5);
   mat B = randu<mat>(4,5);

cx_mat C = zeros<cx_mat>(4,5);

C.set_real(A);
C.set_imag(B);

   mat A = randu<mat>(4,5);
   mat B = randu<mat>(4,5);

cx_mat C = cx_mat(A,B);

mat A;
A.set_size(5,10);

vec q;
q.set_size(100);

mat A = randu<mat>(5,10);

A.shed_row(2);
A.shed_cols(2,4);

mat A = randu<mat>(5,5);
cout << A.size() << endl;

A.clear();
cout << A.empty() << endl;

mat A = zeros<mat>(5,10);

A.submat( 0,1, 2,3 )      = randu<mat>(3,3);
A( span(0,2), span(1,3) ) = randu<mat>(3,3);
A( 0,1, size(3,3) )       = randu<mat>(3,3);

mat B = A.submat( 0,1, 2,3 );
mat C = A( span(0,2), span(1,3) );
mat D = A( 0,1, size(3,3) );

A.col(1)        = randu<mat>(5,1);
A(span::all, 1) = randu<mat>(5,1);

mat X = randu<mat>(5,5);

// get all elements of X that are greater than 0.5
vec q = X.elem( find(X > 0.5) );

// add 123 to all elements of X greater than 0.5
X.elem( find(X > 0.5) ) += 123.0;

// set four specific elements of X to 1
uvec indices;
indices << 2 << 3 << 6 << 8;

X.elem(indices) = ones<vec>(4);

cube A = randu<cube>(2,3,4);

mat B = A.slice(1);

A.slice(0) = randu<mat>(2,3);
A.slice(0)(1,2) = 99.0;

A.subcube(0,0,1,  1,1,2)             = randu<cube>(2,2,2);
A( span(0,1), span(0,1), span(1,2) ) = randu<cube>(2,2,2);
A( 0,0,1, size(2,2,2) )              = randu<cube>(2,2,2);

// add 123 to all elements of A greater than 0.5
A.elem( find(A > 0.5) ) += 123.0;

mat A = zeros<mat>(4,5);
mat B =  ones<mat>(6,7);

A.swap(B);

mat X = randu<mat>(5,5);
X.swap_rows(0,4);

mat A = randu<mat>(4,5);
mat B = A.t();

// C++11 only example

mat A = ones<mat>(4,5);

// add 123 to every element
A.transform( [](double val) { return (val + 123.0); } );

mat A = randu<mat>(5,10);
A.zeros();      // sets all elements to zero
A.zeros(10,20); // sets the size to 10 rows and 20 columns
                // followed by setting all elements to zero

mat A = eye<mat>(5,5);
mat B = 123.0 * eye<mat>(5,5);

vec v = linspace<vec>(10, 20, 5);
mat X = linspace<mat>(10, 20, 5);

vec  v = ones<vec>(10);
uvec u = ones<uvec>(11);
mat  A = ones<mat>(5,6);
cube Q = ones<cube>(5,6,7);

mat  B = 123.0 * ones<mat>(5,6);

imat A = randi<imat>(5, 6);

imat A = randi<imat>(6, 7, distr_param(-10, +20));

arma_rng::set_seed_random();  // set the seed to a random value

vec  v = randu<vec>(5);
mat  A = randu<mat>(5,6);
cube Q = randu<cube>(5,6,7);

arma_rng::set_seed_random();  // set the seed to a random value

mat A = randu<mat>(2, 3);

mat B = repmat(A, 4, 5);

sp_mat A = speye<sp_mat>(5,5);

sp_mat A = sprandu<sp_mat>(100, 200, 0.1);

sp_mat B = spones(A);

sp_mat A = sprandu<sp_mat>(100, 200, 0.1);

vec A = randu<vec>(5);
mat X = toeplitz(A);
mat Y = circ_toeplitz(A);

vec  v = zeros<vec>(10);
uvec u = zeros<uvec>(11);
mat  A = zeros<mat>(5,6);
cube Q = zeros<cube>(5,6,7);

mat A = randu<mat>(5,5);
mat B = abs(A); 

cx_mat X = randu<cx_mat>(5,5);
mat    Y = abs(X);

mat A = randu<mat>(5,5);

mat B = clamp(A, 0.2, 0.8); 

mat C = clamp(A, A.min(), 0.8); 

mat D = clamp(A, 0.2, A.max()); 

cx_mat X = randu<cx_mat>(5,5);
cx_mat Y = conj(X);

mat  A = randu<mat>(5,5);
fmat B = conv_to<fmat>::from(A);

typedef std::vector<double> stdvec;

stdvec x(3);
x[0] = 0.0; x[1] = 1.0;  x[2] = 2.0;

colvec y = conv_to< colvec >::from(x);
stdvec z = conv_to< stdvec >::from(y); 

mat A = randu<mat>(4,5);
mat B = eps(A);

cx_mat C = randu<cx_mat>(5,5);

mat    A = imag(C);
mat    B = real(C);

mat A = randu<mat>(5,5);
mat B = exp(A);

mat X = randu<mat>(5,5);
mat Y = cos(X);

mat A = randu<mat>(5,5);
double x = accu(A);

mat B = randu<mat>(5,5);
double y = accu(A % B);

// operator % performs element-wise multiplication,
// hence accu(A % B) is a "multiply-and-accumulate" operation

rowvec r = randu<rowvec>(5);
colvec q = randu<colvec>(5);
mat    X = randu<mat>(5,5);

// examples of some expressions
// for which optimised implementations exist

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

mat    A = randu<mat>(5,5);
double c = cond(A);

mat    A = randu<mat>(5,5);
double x = det(A);

vec a = randu<vec>(10);
vec b = randu<vec>(10);

double x = dot(a,b);

mat A = randu<mat>(5,5);
mat B = randu<mat>(5,5);

B(1,1) = datum::inf;

cout << is_finite(A) << endl;
cout << is_finite(B) << endl;

cout << is_finite( 0.123456789 ) << endl;
cout << is_finite( datum::nan  ) << endl;
cout << is_finite( datum::inf  ) << endl;

mat A = randu<mat>(5,5);

double val;
double sign;

log_det(val, sign, A);

vec    q = randu<vec>(5);

double x = norm(q, 2);
double y = norm(q, "inf");

mat   A = randu<mat>(4,5);

uword r = rank(A);

mat    A = randu<mat>(5,5);

double x = trace(A);

vec V = randu<vec>(10);
mat X = randu<mat>(5,5);


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

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

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

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

vec V = randu<vec>(10);
mat X = randu<mat>(5,5);


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

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

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

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

mat A = randu<mat>(5,5);

vec d = diagvec(A);

colvec v   = randu<colvec>(10,1);
double val = max(v);

mat    X = randu<mat>(10,10);
rowvec r = max(X);

// same result as max(X)
// the 0 explicitly indicates "traverse across rows"
rowvec s = max(X,0); 

// the 1 explicitly indicates "traverse across columns"
colvec t = max(X,1);

// find the overall maximum value
double y = max(max(X));

// element-wise maximum
mat A = randu<mat>(5,6);
mat B = randu<mat>(5,6);
mat C = arma::max(A,B);  // use arma:: prefix to distinguish from std::max()

colvec q = randu<colvec>(10,1);
double x = prod(q);

mat    A = randu<mat>(10,10);
rowvec b = prod(A);

// same result as prod(A)
// the 0 explicitly indicates
// "traverse across rows"
rowvec c = prod(A,0);

// the 1 explicitly indicates
// "traverse across columns"
colvec d = prod(A,1);

// find the overall product
double y = prod(prod(A));

colvec q = randu<colvec>(10,1);
double x = sum(q);

mat    A = randu<mat>(10,10);
rowvec b = sum(A);

// same result as sum(A)
// the 0 explicitly indicates "traverse across rows"
rowvec c = sum(A,0);

// the 1 explicitly indicates "traverse across columns"
colvec d = sum(A,1);

// find the overall sum
double y = sum(sum(A));

double z = accu(A);

mat A    = randu<mat>(5,5);
mat B    = mean(A);
mat C    = var(A);
double m = mean(mean(A));

vec    q = randu<vec>(5);
double v = var(q);

vec A = randu<vec>(128) - 0.5;
vec B = randu<vec>(128) - 0.5;

vec C = conv(A,B);

mat X = randu<mat>(4,5);
mat Y = randu<mat>(4,5);

mat R = cor(X,Y);
mat S = cor(X,Y, 1);

mat X = randu<mat>(4,5);
mat Y = randu<mat>(4,5);

mat C = cov(X,Y);
mat D = cov(X,Y, 1);

vec a = randu<vec>(3);
vec b = randu<vec>(3);

vec c = cross(a,b);

mat A = randu<mat>(5,5);
mat B = cumsum(A);
mat C = cumsum(A, 1);

vec x = randu<vec>(10);
vec y = cumsum(x);

mat A = randu<mat>(5,5);
mat B = diagmat(A);
mat C = A*diagmat(A);

rowvec q = randu<rowvec>(5);
colvec r = randu<colvec>(5);
mat    X = diagmat(q)*diagmat(r);

mat  A  = randu<mat>(5,5);
mat  B  = randu<mat>(5,5);

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

mat A = randu<mat>(5,5);

A(1,1) = datum::inf;

// accumulate only finite elements
double val = accu( A.elem( find_finite(A) ) );

mat A = randu<mat>(5,5);

A(1,1) = datum::inf;

// change non-finite elements to zero
A.elem( find_nonfinite(A) ).zeros();

mat A = randu<mat>(5,5);

mat B = fliplr(A);
mat C = flipud(A);

vec  v  = randn<vec>(1000); // Gaussian distribution

uvec h1 = hist(v, 11);
uvec h2 = hist(v, linspace<vec>(-2,2,11));

vec  v = randn<vec>(1000);  // Gaussian distribution

uvec h = histc(v, linspace<vec>(-2,2,11));

mat X = randu<mat>(4,5);
mat Y = randu<mat>(20000,30000);

inplace_trans(X);            // use standard algorithm by default

inplace_trans(Y, "lowmem");  // use low-memory (and slow) algorithm

mat A = randu<mat>(4,5);
mat B = randu<mat>(4,6);
mat C = randu<mat>(6,5);

mat X = join_rows(A,B);
mat Y = join_cols(A,C);

mat A = randu<mat>(4,5);
mat B = randu<mat>(5,4);

mat K = kron(A,B);

vec A = randu<vec>(10);
vec B = normalise(A);
vec C = normalise(A, 1);

mat X = randu<mat>(5,6);
mat Y = normalise(X);
mat Z = normalise(X, 2, 1);

mat A = randu<mat>(10, 5);
mat B = reshape(A, 5, 10);

mat A = randu<mat>(4, 5);
mat B = resize(A, 7, 6);

mat A = randu<mat>(4,5);
mat B = shuffle(A);

mat A = randu<mat>(10,10);
mat B = sort(A);

vec  q       = randu<vec>(10);
uvec indices = sort_index(q);

mat A = randu<mat>(5,5);

mat B = symmatu(A);
mat C = symmatl(A);

mat A = randu<mat>(5,10);

mat B = trans(A);
mat C = A.t();    // equivalent to trans(A), but more compact

mat A = randu<mat>(5,5);

mat U = trimatu(A);
mat L = trimatl(A);

mat X;
X << 1 << 2 << endr
  << 2 << 3 << endr;

mat Y = unique(X);

mat A = randu<mat>(4, 5);

vec v = vectorise(A);

mat X = randu<mat>(5,5);
mat Y = X.t()*X;

mat R = chol(Y);

// for matrices with real elements

mat A = randu<mat>(50,50);
mat B = A.t()*A;  // generate a symmetric matrix

vec eigval;
mat eigvec;

eig_sym(eigval, eigvec, B);


// for matrices with complex elements

cx_mat C = randu<cx_mat>(50,50);
cx_mat D = C.t()*C;

   vec eigval2;
cx_mat eigvec2;

eig_sym(eigval2, eigvec2, D);

mat A = randu<mat>(10,10);

cx_vec eigval;
cx_mat eigvec;

eig_gen(eigval, eigvec, A);

mat A = randu<mat>(10,10);
mat B = randu<mat>(10,10);

cx_vec eigval;
cx_mat eigvec;

eig_pair(eigval, eigvec, A, B);

// generate sparse matrix
sp_mat A = sprandu<mat>(1000, 1000, 0.1);
sp_mat B = A.t()*A;

vec eigval;
mat eigvec;

eigs_sym(eigval, eigvec, B, 5);  // find 5 eigenvalues/eigenvectors

// generate sparse matrix
sp_mat A = sprandu<sp_mat>(1000, 1000, 0.1);  

cx_vec eigval;
cx_mat eigvec;

eigs_gen(eigval, eigvec, A, 5);  // find 5 eigenvalues/eigenvectors

   vec X = randu<vec>(100);
cx_vec Y = fft(X, 128);

   mat A = randu<mat>(100,100);
cx_mat B = fft2(A);
cx_mat C = fft2(A, 128, 128);

mat A = randu<mat>(5,5);

mat B = inv(A);

mat A   = randu<mat>(5,5);

mat AtA = A.t() * A;

mat B   = inv_sympd(AtA);

mat A = randu<mat>(5,5);

mat L, U, P;

lu(L, U, P, A);

mat B = P.t()*L*U;

mat A = randu<mat>(4,5);

mat B = pinv(A);        // use default tolerance

mat C = pinv(A, 0.01);  // set tolerance to 0.01