NABLA  Nabla Ain't Basic Linear Algebra
Quick Start Guide

Vectors, which the library is all about, are represented by the vector template class.

In general the usage of vector is very similar to the STL vector or valarray classes. Actually it is closer to valarray which models sequence or C-style array, not the container concept defined by the STL.

Here is how you can create vectors:

vector<double> v;

Here an empty vector is created: it has zero size and holds no elements. Notice that the template must be parameterized with the type of elements, just like the STL vector.


vector<double> v(3);

Here an R3 vector is created (since double models real values). That is a vector of real (double precision floating point) values with size equal to three (3 elements). The size of a given vector can be acquired using size() member function:

const size_t s = v.size();

vector<double> v(3, 0.);

Here a 3 element vector is created with all element values initialized to zero (0.).


To resize a vector to some arbitrary size use the resize() member function:

v.resize(s + 3);

Elements of a vector can be accessed in two ways:

v[0] = 1.;
v(1) = 2.;

The first is the overloaded subscript operator[] and the second is the overloaded function call operator().


Finally a vector can be constructed from as well as assigned a valid vector expression (see documentation for vector_expression ):

vector<double>
sum(v1 + v2),
diff;
diff = v1 - v2;

Note that v1 and v2 must have the same sizes.


There is a set of operations which can be applied to vectors, e.g. sum of two vectors (as above) or difference (as above). There are some other natural operations defined as well, for brief list of them see the Showcase.
Contents of two vector objects can be very efficiently swapped:

swap(v1, v2);
v1.swap(v2);

These two ways have identical effect.


There are also some convenient operations defined. For example a sum of products of corresponding elements of two vectors can be computed:

sum_of_products(v1, v2);

A linear map is a very important term in linear algebra. It is so important that this library dare to implement it in the form of matrices. A matrix is a plain table of values arranged in rows and columns. Formally it can be viewed as a vector of vectors but this is not quite the idea. In this library it is represented by the matrix template class.

To create a rectangular matrix write something like

matrix<double> m;

Here m is an empty matrix with no elements. As in the vector case you can specify the size of the matrix and a value for initialization:

matrix<double>
m1(3, 3), //specifying size
m2(2, 3, 1.); //specifying size and value

The number of rows and columns of a matrix can be acquired using rows() and cols() member functions:

const size_t
rows = m.rows(),
columns = m.cols();

A matrix can be resized using resize() member function:

m.resize(rows + 1, columns);

The most efficient way (in run time view) to access an element is

m(0, 1) = 1.;

using overloaded function call operator().

However there are numerous alternative ways to do that:

m[0][0] = 1.; //classics
m[0](1) = 2.;
m(1)[0] = 3.;
m(1)(1) = 4.;

Pick your favourite!


A matrix can be constructed from and assigned a valid matrix expression (see documentation for matrix_expression ):

matrix<double>
sum(m1 + m2),
prod;
prod = m1*m2;

Note that m1 and m2 must have conforming sizes.


There is a set of operations which can be applied to matrices, e.g. sum of two matrices (as above) or product (as above). There are some other natural operations defined as well, for brief list of them see the Showcase.
Contents of two matrix objects can be very efficiently swapped:

swap(m1, m2);
m1.swap(m2);

These two ways have identical effect.


There are also some convenient operations defined. For example rows of a matrix can be swapped:

swap(row(m, 0), row(m, 1));

A matrix is a linear map of one Rn vector to another Rm vector. The map operation is simply computing weighted sum of the source vector elements using matrix rows as weighting coefficients. Thus a matrix is said to be a m rows and n columns matrix, or simply an m by n matrix. If m is equal to n then such a matrix is said to be square. In pseudocode the map operation performs like this:

for each i
{ vm[i] = 0;
for each j
vm[i] += m[i][j]*v[j]; }

Here v is the source vector, vm is the mapped vector and m is the mapping matrix.

To map a vector using a matrix use the map() function:

vm = map(m, v);

It is very convenient to use matrix algebra to perform map operations. This library provides it the natural way:

vm = m*v;

Here m*v represents a well defined matrix-vector product.


Summing it all up we can write a simple (naive) matrix inversion algorithm:

template<typename value_t>
{//in-place matrix inversion algorithm using Gauss-Jordan elimination
const size_t n = m.rows();
if (n==0 || n!=m.cols())
throw "matrix has wrong size";
for (size_t i = 0; i<n; ++i) //iterating over columns
{
const value_t pivot = m(i, i); //so called pivot element
if (pivot==traits<value_t>::zero)
throw "matrix is singular";
m(i, i) = traits<value_t>::unit; //substituting identity matrix into 'm'
row(m, i) *= value_t(1)/pivot; //dividing row with pivot element
for (size_t j = 0; j<n; ++j) //iterating over all other rows
{
if (i==j) //we don't need row with pivot element
continue;
const value_t e = m(j, i); //multiplier for the row with pivot
m(j, i) = traits<value_t>::zero; //substituting identity matrix into 'm'
row(m, j) -= e*row(m, i); //eliminating row
} //^^the third cycle is actually here so this is an O(n^3) operation
}
return m;
}

Personally I consider this implementation both concise and conscious. The writer's intention is explicitly expressed and all mess is hidden away.