So, I'm testing the Armadillo matrix library that undergirds mlpack. I'm certainly no math guy, but this example code they provide in the archive seems to be doing a lot of work. I added simple timer code and it appears to be coming in consistently at 40-45milliseconds to complete. That's with all the stream outputs to the console.
#include <chrono>
#include <iostream>
#include <armadillo>
using namespace std;
using namespace arma;
using chrono::milliseconds;
using chrono::steady_clock;
int main(int argc, char** argv) {
cout << "Armadillo version: " << arma_version::as_string() << '\n';
steady_clock clock{};
auto begin = clock.now();
mat A(2, 3); // directly specify the matrix size (elements are uninitialised)
// .n_rows and .n_cols are read only
cout << "A.n_rows: " << A.n_rows << '\n';
cout << "A.n_cols: " << A.n_cols << '\n';
A(1, 2) = 456.0; // directly access an element (indexing starts at 0)
A.print("A:");
A = 5.0; // scalars are treated as a 1x1 matrix
A.print("A:");
A.set_size(4, 5); // change the size (data is not preserved)
A.fill(5.0); // set all elements to a particular value
A.print("A:");
A = {{0.165300, 0.454037, 0.995795, 0.124098, 0.047084},
{0.688782, 0.036549, 0.552848, 0.937664, 0.866401},
{0.348740, 0.479388, 0.506228, 0.145673, 0.491547},
{0.148678, 0.682258, 0.571154, 0.874724, 0.444632},
{0.245726, 0.595218, 0.409327, 0.367827, 0.385736}};
A.print("A:");
// determinant
cout << "det(A): " << det(A) << '\n';
// inverse
cout << "inv(A): " << '\n' << inv(A) << '\n';
// save matrix as a text file
A.save("A.txt", raw_ascii);
// load from file
mat B;
B.load("A.txt");
// submatrices
cout << "B( span(0,2), span(3,4) ):" << '\n'
<< B(span(0, 2), span(3, 4)) << '\n';
cout << "B( 0,3, size(3,2) ):" << '\n' << B(0, 3, size(3, 2)) << '\n';
cout << "B.row(0): " << '\n' << B.row(0) << '\n';
cout << "B.col(1): " << '\n' << B.col(1) << '\n';
// transpose
cout << "B.t(): " << '\n' << B.t() << '\n';
// maximum from each column (traverse along rows)
cout << "max(B): " << '\n' << max(B) << '\n';
// maximum from each row (traverse along columns)
cout << "max(B,1): " << '\n' << max(B, 1) << '\n';
// maximum value in B
cout << "max(max(B)) = " << max(max(B)) << '\n';
// sum of each column (traverse along rows)
cout << "sum(B): " << '\n' << sum(B) << '\n';
// sum of each row (traverse along columns)
cout << "sum(B,1) =" << '\n' << sum(B, 1) << '\n';
// sum of all elements
cout << "accu(B): " << accu(B) << '\n';
// trace = sum along diagonal
cout << "trace(B): " << trace(B) << '\n';
// generate the identity matrix
mat C = eye<mat>(4, 4);
// random matrix with values uniformly distributed in the [0,1] interval
mat D = randu<mat>(4, 4);
D.print("D:");
// row vectors are treated like a matrix with one row
rowvec r = {0.59119, 0.77321, 0.60275, 0.35887, 0.51683};
r.print("r:");
// column vectors are treated like a matrix with one column
vec q = {0.14333, 0.59478, 0.14481, 0.58558, 0.60809};
q.print("q:");
// convert matrix to vector; data in matrices is stored column-by-column
vec v = vectorise(A);
v.print("v:");
// dot or inner product
cout << "as_scalar(r*q): " << as_scalar(r * q) << '\n';
// outer product
cout << "q*r: " << '\n' << q * r << '\n';
// multiply-and-accumulate operation (no temporary matrices are created)
cout << "accu(A % B) = " << accu(A % B) << '\n';
// example of a compound operation
B += 2.0 * A.t();
B.print("B:");
// imat specifies an integer matrix
imat AA = {{1, 2, 3}, {4, 5, 6}, {7, 8, 9}};
imat BB = {{3, 2, 1}, {6, 5, 4}, {9, 8, 7}};
// comparison of matrices (element-wise); output of a relational operator is a
// umat
umat ZZ = (AA >= BB);
ZZ.print("ZZ:");
// cubes ("3D matrices")
cube Q(B.n_rows, B.n_cols, 2);
Q.slice(0) = B;
Q.slice(1) = 2.0 * B;
Q.print("Q:");
// 2D field of matrices; 3D fields are also supported
field<mat> F(4, 3);
for (uword col = 0; col < F.n_cols; ++col)
for (uword row = 0; row < F.n_rows; ++row) {
F(row, col) = randu<mat>(2, 3); // each element in field<mat> is a matrix
}
F.print("F:");
auto end = clock.now();
cout << duration_cast<milliseconds>(end - begin).count() << "ms\n";
return 0;
}
Maybe one of you math wizards can look this over and see if you think it's doing well for yourselves. Will post outputs next.