• Description
• Classes
• Element Access
• Arithmetic Operators and Functions
• Constants
• Constructors
• Input and Output
• Transformations
• Sparse Vectors and Matrices
• Sub Vectors
• Solving Systems of Linear Equations
• Factoring Matrices
• Compiling with VL
• Using VL With OpenGL
• Author

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;

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]

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

    (*vPtr)[20] = 3.0;

    Float *vecDataPtr = v.Ref(), *matDataPtr = m.Ref();

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

Arithmetic Operators and Functions

Arithmetic Operators

    Basic arithmetic:         + - * /
    Accumulation arithmetic:  += -= *= /=
    Comparison:               ==, !=

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

    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]

Matrix Functions

    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.

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]

    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

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

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

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

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;

    [1 1 1][2 2 2]

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

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

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

    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]

Manipulation

*    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

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

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

    sv1.Set(10, sv2.Get(20)); // sv1[10] = sv2[20]

Sparse Fuzziness

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.

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.

    diag(m) *= 2.0;
    sub(m, 2, 2) += Matf(2, 2, vl_1);

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 + A^T)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

The solvers return the squared residual of the linear system, ||Ax - b||², 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);

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 = UDV^T, where:

• U is m x n, and orthogonal.
• D is n x n, and is diagonal and positive semi-definite; its elements are the square roots of the eigenvalues of A^TA.
• V is n x n, and orthogonal.

• The SVD says that we can view any matrix transformation as: rotate the vector from source space, scale about the axes, and rotate the vector again into the destination space.
•  The matrix's condition number is ratio of the largest entry of D to the smallest.
• U and V^T are column-orthonormal.
•  The smallest least squares solution of Ax = b is x = V F U^T b, where F = D^-1, but with zero entries on the diagonal of D kept as zero rather than replaced with their inverse.

    #include "Factor.h"

    Void QRFactorization(Matd &A, Matd &Q, Matd &R);
    Void SVDFactorization(Matd &A, Matd &U, Matd &V, Vecd &diagonal);

Headers

Linking

Compile options

    VL_CHECKING     - turn on index checking and assertions
    VL_ROW_ORIENT   - transformations operate on row vectors instead of column vectors

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

Please forward bug reports, comments, or suggestions to: