VL: Vector Library        (1.2)
Contents
Description

The Vector Library (VL) provides a set of vector and matrix classes, as well as a number of functions for performing arithmetic with them. Equation-like syntax is supported via C++ class operators, for example:

   
    #include "VLfd.h"
    Vec3f   v(1.0, 2.0, 3.0);
    Mat3d   m(1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0);

    v =  2 * v + m * v;
    v *= (m / 3.0) * norm(v);
    
    cout << v << endl;
Both generic (arbitrarily-sized), and fixed-size (2, 3 and 4 element) vectors and matrices are supported. The latter are provided for the efficient manipulation of vectors or points in 2D or 3D space, and make heavy use of inlining for frequently-used operations. (One of the design goals of VL was to ensure that it was as fast as the C-language, macro-based libraries it was written to replace.)

Vectors and matrices can be composed of either floats or doubles; the element type is indicated by the suffix. It is possible to mix (for example) matrices of doubles with vectors of floats, as in the example above. It is also possible to instantiate VL for other other element types with their own suffixes (e.g., complex numbers).

VL also contains classes for sparse vector/matrices, sub-vector/matrices, and implementations of some iterative solvers.

VL is free for commercial and non-commercial use; see the LICENSE file for redistribution conditions. (These apply only to the source code; binaries may be freely redistributed, no strings attached.)

VL requires C++. It is known to compile under g++, MSVC++ 5.0, Irix CC, and Metrowerks C++ (macintosh). The latest version can be retrieved from http://www.cs.cmu.edu/~ajw/public/dist/. This documentation can be found online at http://www.cs.cmu.edu/~ajw/doc/vl.html.
 
 
Classes

VL contains the following types and classes:

    Fixed-size:
        Vec2[fd]        2-vector
        Vec3[fd]        3-vector
        Vec4[fd]        4-vector
        Mat2[fd]        2 x 2 matrix
        Mat3[fd]        3 x 3 matrix
        Mat4[fd]        4 x 4 matrix
    Generic:
        Vec[fd]         n-vector
        Mat[fd]         n x m matrix
    Sparse:
        SparseVec[fd]   n-vector optimised for sparse storage
        SparseMat[fd]   n x m matrix optimised for sparse storage
    Sub:
        SubVec[fd]      n-vector which is a subset of another vector
        SubMat[fd]      n x m matrix which is a subset of another matrix
        SubSVec[fd]     the same for sparse vectors & matrices
        SubSMat[fd]
 
Element Access

The elements of a vector or matrix are accessed with standard C array notation:

   
    v[2] = 4.0;         // set element 2 of the vector
    m[3][4] = 5.0       // set row 3, column 4 of the matrix
    m[2] = v;           // set row 2 of the matrix
For the resizeable vector types, the current size can be obtained from the Elts() method for vectors, and the Rows() and Cols() methods for matrices. To iterate over all members of these types, you can use code of the form:
    for (i = 0; i < v.Elts(); i++)
        v[i] = i;
    for (i = 0; i < m.Rows(); i++)
        for (j = 0; j < m.Cols(); j++)
            m[i][j] = i + j;
Though it seems slightly unintuitive, if you have a pointer to a vector or matrix, you must dereference it first before indexing it:
    (*vPtr)[20] = 3.0;
If you need a pointer to the data belonging to a vector or matrix, use the Ref() method. (Matrices are stored by row.)
    Float *vecDataPtr = v.Ref(), *matDataPtr = m.Ref();
Warning: Any pointer to the data of a generic matrix or vector will become invalid as soon as it is resized.

Note: If you compile with -DVL_CHECKING, index range checks will be performed on all element accesses.
 
 
Arithmetic Operators and Functions

Arithmetic Operators

The following binary operators are defined for all vector and matrix classes, as long as both operands are of the same type.
    Basic arithmetic:         + - * /
    Accumulation arithmetic:  += -= *= /=
    Comparison:               ==, !=
Vector multiplication and division is pairwise: (a * b)[i] = a[i] * b[i]. (See below for how to form the dot product of two vectors with dot().) Matrix multiplication is defined as usual, and matrix division is undefined.

For both matrices and vectors, multiplication and division by a scalar is also allowed. Matrices can be multiplied either on the left or the right by a vector. In the expression m * v, v is treated as a column vector; in the expression v * m, it is treated as a row vector.

Vector Functions

The following is a list of the various vector functions, together with a short description of what they return.
    Float   dot(const Vec[fd] &a, const Vecf &b);  // inner product of a and b
    Float   len(const Vecf &v);                    // length of v: || v ||
    Float   sqrlen(const VecNf &v);                // length of v, squared
    VecNf   norm(const VecNf &v);                  // v / || v ||

    Vec2f   cross(const Vec2f &a);                 // vector orthogonal to a 
    Vec3f   cross(const Vec3f &a, const Vec3f &b); // vector orthogonal to a and b
    Vec4f   cross(const Vec4f &a, const Vec4f &b, const Vec4f &c);
                                                   // vector orthogonal to a, b and c

    Vec2f   proj(const Vec3f &v);                  // homog. projection: v[0..1] / v[2]
    Vec3f   proj(const Vec4f &v);                  // homog. projection: v[0..2] / v[3]
In the above, VecN is either a Vec or a Vec[234], and all functions have corresponding Double/VecNd versions. For more on the use of the proj() operator, see Transformations.

Matrix Functions

The following functions can be used with matrices.
    MatNf   trans(const MatNf &m);                  // Transpose of m           
    Float   trace(const MatNf &m);                  // Trace of m
    MatNf   adj(const MatNf &m);                    // Adjoint of m
    Float   det(const MatNf &m);                    // Determinant of m
    MatNf   inv(const MatNf &m);                    // Inverse of m, if it exists.
Here MatN is any matrix type, though the det() and adj() functions are only defined for Mat[234][fd].
 
 
Constants

There are a number of 'magic' constants in VL that can be used to initialise vectors or matrices with simple assignment statements. For example:

    Vec3f v; Mat3f m; Vecf v8(8);

    v = vl_0            [0, 0, 0]
    v = vl_y            [0, 1, 0]
    v = vl_1            [1, 1, 1]
    
    m = vl_0;           3 x 3 matrix, all elts. set to zero.
    m = vl_1;           3 x 3 identity matrix
    m = vl_B;           3 x 3 matrix, all elts. set to one.
    
    v8 = vl_axis(6);    [0, 0, 0, 0, 0, 0, 1, 0]
Below is a summary of the constants defined by VL.
    vl_one/vl_1/vl_I              vector of all 1s, or identity matrix
    vl_zero/vl_0/vl_Z             vector or matrix of all 0s
    vl_B                          matrix of all 1s
    vl_x, vl_y, vl_z, vl_w        x, y, z and w axis vectors
    vl_axis(n)                    zero vector with element n set to 1
    vl_pi                         pi!
    vl_halfPi                     pi/2
 
Constructors

In general, a vector or matrix constructor should be given either one of the initialiser constants listed above, or a list of values for its elements. If neither of these is supplied, the variable will be uninitialised. The first arguments to the constructor of a generic vector or matrix should always be the required size. Thus matrices and vectors are declared as follows:

    Vec[fd][234]    v([initialisation_constant | element_list]);
    Vec[fd]         v([elements, [initialisation_constant | element_list]]);
    Mat[fd][234]    m([initialisation_constant | element_list]);
    Mat[fd]         m([rows, columns, [initialisation_constant | element_list]]);
If generic vectors or matrices are not given a size when first created, they are regarded as empty, with no associated storage. This state persists until they are assigned a matrix/vector or the result of some computation, at which point they take on the dimensions of that result.

Examples:

    Vec3f    v(vl_1);
    Vec3f    v(1.0, 2.0, 3.0);
    Vecf     v(6, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0); 
    Vecf     v(20, vl_axis(10));
    Mat2f    m(1.0, 2.0, 3.0, 4.0);
    Matf     m(10, 20, vl_I);
Warning: When initialising a generic vector or matrix with a list of elements, you must always ensure there is no possibility of the element being mistaken for an integer. (This is due to limitations of the stdarg package.) Make sure that each element value has either an exponent or a decimal point, i.e., use '2.0' rather than just '2'.

Finally, to set the size of a empty matrix or vector explicitly, or resize an existing matrix or vector, use the SetSize method:

    v.SetSize(23);
    m.SetSize(10, 20);
 
Input and Output

All of the vector and matrix types in VL can be used in iostream-type expressions. For example:

    #include <iostream.h>
    Vec3d v(vl_1);
    cout << v << 2 * v << endl;
    cin >> v;
will output
    [1 1 1][2 2 2]
and then prompt for input. Vectors and matrices are parsed in the same format that they are output: vectors are delimited by square brackets, elements separated by white space, and matrices consist of a series of row vectors, again delimited by square brackets.
 
 
Transformations

The following are the transformations supported by VL.

    Mat2f   Rot2f(Double theta)
            // rotate a 2d vector CCW by theta
    Mat2f   Scale2f(const Vec2f &s)
            // scale by s around the origin
    
    Mat3f   HRot3f(Double theta)
            // rotate a homogeneous 2d vector CCW by theta
    Mat3f   HScale3f(const Vec2f &s)
            // scale by s around the origin, in homogeneous 2d coords.
    Mat3f   HTrans3f(const Vec2f &t)
            // translate a homogeneous 2d vector by t

    Mat3f   Rot3f(const Vec3f &axis, Double theta)
            // rotate a 3d vector CCW around axis by theta
    Mat3f   Rot3f(const Vec4f &q)
            // rotate a 3d vector by the quaternion q
    Mat3f   Scale3f(const Vec3f &s)
            // scale by s around the origin

    Mat4f   HRot4f(const Vec3f &axis, Double theta)
            // rotate a homogeneous 3d vector CCW around axis by theta
    Mat4f   HRot4f(const Vec4f &q)
            // rotate a homogeneous 3d vector by the quaternion q
    Mat4f   HScale4f(const Vec3f &s)
            // scale by s around the origin, in homogeneous 3d coords
    Mat4f   HTrans4f(const Vec3f &t)
            // translate a homogeneous 3d vector by t
All transformations have corresponding Double versions with a 'd' suffix. Transformations with a prefix of 'H' operate in the homogeneous coordinate system, which allows translation and shear transformations, as well as the usual rotation and scale. In this coordinate system an n-vector is embedded in a (n+1)-dimensional space, e.g., a homogeneous point in 2d is represented by a 3-vector.

To convert from non-homogeneous to homogeneous vectors, make the extra coordinate (usually 1) the second argument in a constructor of/cast to the next-higher dimension vector. To project from a homogeneous vector down to a non-homogeneous one (performing a homogeneous divide in the process), use the proj() function. This process can be simplified by the use of the xform() function, which applies a transform to a vector or composes two transformations, performing homogeneous/nonhomogeneous conversions as necessary. For example:

    Vec3d x,y;

    // apply homogeneous transformations to a 3-vector, assuming column vectors 
    x = proj(Scale4d(...) * Trans4d(...) * Vec4d(y, 1.0));
    // do the same thing with xform()
    x = xform(xform(Scale4d(...), Trans4d(...)), y);
By default, VL assumes that transformations should operate on column vectors (v = T * v), though it can be compiled to assume row vectors instead (v = v * T). Using the xform functions is a good way of isolating yourself from this assumption.
 
 
Sparse Vectors and Matrices

VL contains both a sparse vector type, which stores only the non-zero elements of the vector, and a sparse matrix type, whose rows are sparse vectors. Sparse vectors can be efficiently combined with other sparse vectors or normal vectors using the standard vector operations:

    SparseMatf sm;
    SparseVecf sv1, sv2, sv3;
    Vecf  v;
    v = sm * v;
    sm[0] += v;
    sv1 = sv2 + sv3;

Initialisation

Sparse vectors are typically initialised by giving a list of index, element pairs:
    SparseVecf sv;      // unsized sparse vector
    SparseVecf sv(20);  // zero vector of length 20.
    SparseVecf sv(5, 1, 1.0, 4, 4.0, VL_SV_END);
                        // the vector [0.0, 1.0, 0.0, 0.0 4.0]
The standard vector initialisers can also be used.
 

Manipulation

Once you have your sparse vector, it can be changed in the following ways:

*    Re-initialise with the SetElts() method:

    sv.SetSize(5);
    sv.SetElts(1, 1.0, 4, 4.0, VL_SV_END);
    // sets sv to [0.0, 1.0, 0.0, 0.0 4.0]


*    Use the SVIter[fd] iterator, which lets you iterate over the elements of a sparse vector, and access them using the methods:

j.Data() : returns the current element's data
j.Index() : returns the current element's index
Note:  It is highly recommended that you use the SVIter[fd] class to manipulate sparse vectors, as it is written to be as efficient as possible, even performing a binary search for elements when necessary. The iterator class can be used in a number of ways:
* Use Begin(), Inc(), AtEnd() to iterate over the non-zero elements of the vector:
    SVIterf j;
    // sv = sv * 2
    for (j.Begin(sv); !j.AtEnd(); j.Inc())
        j.Data() *= 2.0;
*   Use one of the following methods:
  Inc(Int i)    moves on to elt i, where i will increase by 1 on each call
  IncTo(Int i)  moves on to elt i, where i will increase monotonically
within another for-loop to access the elements of the sparse vector corresponding to i. For example:
    // v = v + sv
    for (j.Begin(sv), i = 0; i < v.NumItems(); i++)
    {
        j.Inc(i);
        v[i] += j.Data();
    }

    // a += dot(sv1, sv2)
    for (j.Begin(sv2), k.Begin(sv1); !k.AtEnd(); k.Inc())
    {
        j.IncTo(k.Index());  // find corresponding elt in sv2
        if (j.Exists())
            a += j.Data() * k.Data();
    }
*   Use the Overlay() method: a.Overlay(b) performs a[i] = b[i] for all non-zero b[i].

*    Direct access: Begin(),AddElt() or AddNZElt() new element pairs in order, then call End(). (Use AddNZElt() if you know for certain the added element will be non-zero.) For example:

    // set sv to [0.0, 1.0, 0.0, 0.0 4.0]
    sv.NumElts(5);
    sv.Begin();
    sv.AddNZElt(1, 1.0);
    sv.AddNZElt(4, 4.0);
    sv.End();
*   As a last resort, use the Get() and Set() methods. These calls are not efficient for multiple accesses, but will at least perform a binary search to locate the requested element quickly.
    sv1.Set(10, sv2.Get(20)); // sv1[10] = sv2[20]
Note:  The best way to write code for sparse vectors and matrices is to use the SVIter[fd] class, and recast code to use the efficient vector operations where possible.
 

Sparse Fuzziness

The SparseVec class has a tolerance level for elements to be considered zero, referred to as the fuzz. (If |x| < fuzz, it is treated as zero.) This can be set with the method SparseVec[fd]::SetFuzz(fuzz). The default value of fuzz is 1e-10.

A convenient way to test if an element is zero according to the current fuzz setting is to use the SparseVec[fd]::IsNonZero(elt) method.
 
 
Sub Vectors

VL provides the following functions for accessing subregions of vectors and matrices:

    Vec[fd]     sub(const Vec[fd] &v, Int start, Int length);  
    Vec[fd]     first(const Vec[fd] &v, Int length);  
    Vec[fd]     last(const Vec[fd] &v, Int length);
    SubMat[fd]  sub(const Mat[fd] &m, Int top, Int left, Int height, Int width = 1);
    SubMat[fd]  sub(const Mat[fd] &m, Int rows, Int cols);
    SubVec[fd]  col(const Mat[fd] &m, Int i);
    SubVec[fd]  row(const Mat[fd] &m, Int i);
    SubVec[fd]  diag(const Mat[fd] &m, Int diagNum);


The utility of these functions is best illustrated with some examples:

    u = sub(v, 2, 4);         // return the 4 elements of v starting at element 2.
    u = first(v, 2);          // return the first 2 elements of v.
    u = last(v, 2);           // return the last 2 elements of v.

    v = m[i];                 // extract row i of m
    v = col(m, i);            // extract column i of m
    v = row(m, i);            // extract row i of m

    v = diag(m);              // extract main diagonal of m
    v = diag(m, i)            // extract diagonal starting on column i
    v = diag(m, -i)           // extract diagonal starting on row i

    n = sub(m, 2, 3);         // returns the upper-left 2 rows and three columns of m.
    n = sub(m, i, j, 2, 3);   // returns the 2 rows and 3 columns of m starting at
                              // row i, column j.
Warning: remember that indexing is 0-based in VL, so row 2 above refers to the third row from the top of the matrix, and so on.

The subvector and submatrix types returned by the sub(), col() and diag() functions can, in addition to being assigned to normal vectors as above, also be assigned to:
 

    diag(m) = diag(m) * 2.0;     // multiply diagonal elements of m by 2
    sub(m, 2, 2) = sub(m, 2, 2) + Matf(2, 2, vl_1);
                                 // add 1 to each of the upper-left 2 x 2
                                 // elements of m.
Warning: The standard in-place operations are not defined on submatrix regions, so the following will not work:
    diag(m) *= 2.0;
    sub(m, 2, 2) += Matf(2, 2, vl_1);
 
Solving Systems of Linear Equations

VL comes with a number of solvers, which are routines for finding the solution of the linear system of equations Ax = b. Currently these include SolveOverRelax(), which uses the overrelaxation form of Gauss-Seidel iteration, and SolveConjGrad(), which uses the conjugate gradient algorithm. Conjugate gradient is asymptotically faster than Gauss-Seidel, but it assumes that A is both positive definite and symmetric. If A is not symmetric, the routine instead solves the system 0.5(A + AT)x = b.

The solvers are defined as follows:

    Double SolveOverRelax(const [Sparse]Matd &A, Vec[fd] &x, const Vec[fd] &b,
                        Double epsilon, Double omega = 1.0, Int *steps = 0);
    Double SolveConjGrad(const [Sparse]Matd &A, Vec[fd] &x, const Vec[fd] &b,
                        Double epsilon, Int *steps = 0);

Parameters

Each iteration of a solver modifies the current approximate solution x. You must set x to an initial guess before first calling the solver routine; a good starting value is often just b.

The solvers return the squared residual of the linear system, ||Ax - b||2, which is a measure of the error in the solution.

The epsilon parameter controls the accuracy of the solution: the solver will return as soon as its estimate of the squared residual drops below epsilon.

For SolveOverRelax, omega controls the amount of overrelaxation. A value of one gives pure Gauss-Seidel iteration. Values higher than that cause the solver to overshoot towards the estimated solution on each iteration. If the system is smooth and well behaved, this can lead to faster convergence times. Generally, setting omega somewhere between 1 and 2 results in the fastest convergence rate, but the exact value is system-specific.

If you want, you can supply a maximum number of iterations to perform via steps. In this case, the routines will set the actual number of interations performed when they return. This can be useful if you wish to interleave steps of the solver with some other activity.

Examples

    // solve Ax = b from an initial guess of x = b
    x = b;
    SolveOverRelax(A, x, b, 1e-6);
    // perform one iteration of the conjugate gradient solver
    Int steps = 1;
    error = SolveConjGrad(A, x, b, 0, &steps);
 
Factoring Matrices

VL contains two routines for factoring matrices; the QR-factorization, and the SVD (singular value decomposition). The former factors any given matrix A into two matrices Q and R, such that A = QR. The Q matrix is orthogonal, and the R matrix is upper-triangular.

The SVD decomposes an m x n matrix A into three matrices, A = UDVT, where:

The SVD has the following interesting properties: The factoring routines are defined as follows:
    #include "Factor.h"
    Void QRFactorization(Matd &A, Matd &Q, Matd &R);
    Void SVDFactorization(Matd &A, Matd &U, Matd &V, Vecd &diagonal);
Both routines destroy the input matrix, A. Currently, it is required that A have the same or more rows than columns. If your matrix has more columns than rows, add enough zero rows to the bottom of it to make it square.
 
Compiling with VL

Headers

For basic use, the only header file needed is VL.h.

Linking

For your final build, link with -lvl (libvl.a). To use the debugging version of VL, which has assertions and range checking turned on, use -lvl.dbg (libvl.dbg.a), and add -DVL_CHECKING to your compile flags. This debugging version includes checks for correct matrix and vector sizes during arithmetic operations.

Compile options

VL uses the following compile-time options:
    VL_CHECKING     - turn on index checking and assertions
    VL_ROW_ORIENT   - transformations operate on row vectors instead of column vectors
 
Using VL With OpenGL

VL comes with a header file, VLgl.h, which makes using VL vectors with OpenGL more convenient. For example:

    #include "VLgl.h"
    
    Vec3f x(24, 0, 100), y(40, 20, 10);
    
    glBegin(GL_LINES);
    glVertex(x); 
    glVertex(y);
    glEnd();
 
Author

Please forward bug reports, comments, or suggestions to:

Andrew Willmott (ajw+vl@cs.cmu.edu), Graphics Group, SCS, CMU.