Matrix Classes

Table Of Contents

• Fixed-Size Matrices
• Dense Higher-Dimensional Matrices
• Sparse Higher-Dimensional Matrices

Fixed-Size Matrices

The fixed-size matrices come in many flavors: The HDK includes classes for 2x2, 3x3, and 4x4 matrices, with separate versions for single-precision and double-precision floating point numbers. The single-precision matrix classes are UT_Matrix2, UT_Matrix3, and UT_Matrix4. The double-precision versions are UT_DMatrix2, UT_DMatrix3, and UT_DMatrix4. The matrix classes interact with the fixed-size tuple classes UT_Vector2, UT_Vector3, and UT_Vector4.

The discussion below focuses on UT_Matrix4 and its interaction with UT_Vector3 and UT_Vector4. The classes UT_Matrix2 and UT_Matrix3 are mostly special cases of UT_Matrix4. The double-precision UT_DMatrix* versions have the same interface as their single-precision UT_Matrix* counterparts.

The matrix class UT_Matrix4 is typically used to represent a change of coordinates, for example between object space and world space. UT_Matrix4 is also commonly used to represent transformations of objects within the same space, for example, scaling, shearing, translation and rotation.

The following conventions apply to the HDK matrix classes:

• The matrix elements are stored in row-major format.
• Positions and vectors are interpreted as row vectors.

Positions and vectors can both be represented by UT_Vector3. Positions are using the rowVecMult function:

UT_Matrix4 object_to_world;
UT_Vector3 object_position;
...
UT_Vector3 world_position = rowVecMult(object_position, object_to_world);

Please note that the left argument is the position and the right argument is the matrix. When your UT_Vector3 stores a vector (as opposed to a position), then you should use the method rowVecMult3 to transform it:

UT_Matrix4 object_to_world;
UT_Vector3 object_vector;
...
UT_Vector3 world_vector = rowVectMult3(object_vector, object_to_world);

The rowVecMult3 method ignores everything but the upper left 3x3 submatrix, so that the translation part of the matrix is not applied.

It is important to realize that the fourth column of the matrix is ignored (considered zero) when UT_matrix4 operates on UT_Vector3 objects. This is fine when you only need to represent affine transformations such as rotation, scaling, shearing, and translation. However, if you do need to represent perspective transformations, then you should use UT_Matrix4 in combination with UT_Vector4. In that case, UT_Vector4 behaves as a tuple of homogeneous coordinates. An overloaded version of the rowVecMult method can then be used to transform both positions and vectors:

UT_Matrix object_to_world;
...
UT_Vector4 object_position(1, 2, 3, 1); // position (1, 2, 3)
UT_Vector4 object_vector(4, 5, 6, 0); // vector (4, 5, 6)

UT_Vector4 world_position = rowVecMult(object_position, object_to_world);
UT_Vector4 world_vector = rowVecMult(object_vector, object_to_world);

Setting the last coordinate of UT_Vector4 to 1 (the default) will create a point, setting it to 0 will create a vector.

The conventions for the HDK matrix classes imply that matrices work from left to right when they are multiplied:

UT_Matrix4 object_to_world;
UT_Matrix4 world_to_view;
...
UT_Matrix4 object_to_view = object_to_world * world_to_view;

The translate, scale and rotate methods can be used to construct transformation matrices. For example, the following code snippet constructs a transformation that first translates by (1, 2, 3) and then scales by (4, 5, 6):

UT_Matrix4 t(1.0) // the identity
t.translate(1, 2, 3); // t = t * "translation by (1, 2, 3)"
t.scale(4, 5, 6); // t = t * "scaling by (4, 5, 6)"

Calling these members is the same as right-multiplying t by a translation, scale, or rotation matrix. The methods pretranslate, prescale and prerotate are similar, but they have the effect of left-multiplying t.

Transformation matrices can be constructed more explicitly by setting all of the matrix elements. Matrix indices are zero-based. The following code fragment constructs the transformation that scales a point by 0.5 and subsequently translates the point by (1, 2, 3).

UT_Matrix t(1f); // create identity matrix
// set scaling elements
t(0, 0) = 0.5f;
t(1, 1) = 0.5f;
t(2, 2) = 0.5f;
// set translation
t(0, 3) = 1.0f; 
t(1, 3) = 2.0f; 
t(2, 3) = 3.0f; 

The invert method can be used to find the inverse of a matrix without modifying the original matrix:

UT_Matrix4 object_to_world;
...
UT_Matrix4 world_to_object;
if( object_to_world.invert(world_to_object) == 0 )
    // world_to_object is now the inverse of object_to_world
else
    // Couldn't invert: object_to_world was singular

Note that the return value is 0 if invert was successful.

Dense Higher-Dimensional Matrices

UT_Matrix can represent non-square matrices of higher dimensions. The representation is dense; all elements are stored, even if they're zero. UT_Matrix is best used when the the number of rows and columns is small, but higher than 4. Six by six UT_Matrix objects are commonly used in the physics simulation code.

Indexing UT_Matrix can be either 0-based or 1-based. This is controlled by the constructor parameters. The constructor expects inclusive ranges for row numbers and the column numbers. For example, the following code creates a UT_Matrix that has m rows and n columns, using 1-based indices:

UT_Matrix A(1, m, 1, n);

UT_Matrix objects can be multiplied with each other using the * operator. There are also various ways to multiply UT_Matrix object and UT_Vector objects:

UT_Matrix A(1, m, 1, n);
UT_Vector input_row(1, m);
UT_Vector input_column(1, n);
...
UT_Vector result_row(1, n);
// result_row = input_row * A;
A.preMult(input_row, result_row);
...
UT_Vector result_column(1, n);
// result_column = A * input_column;
A.postMult(input_column, result_column);

The UT_MatrixSolver class declares utility methods that solve linear systems of equations that are represented using UT_Matrix. In particular, the LU-decomposition method LUDecomp can be used to express a matrix A as a product P L U of a permutation P, a unit lower-triangular matrix L and an upper-triangular matrix U. Once this decomposition has been computed, it can be re-used to efficiently solve multiple linear systems A x = b that involve the same A.

The LUDecomp method stores the resulting matrices L and U in a single UT_Matrix. This is possible because the unit elements on the diagonal of L don't need to be explicitly stored. The associated permutation matrix P is stored in a UT_Permutation object, which takes up much less space than a UT_Matrix. The method LUBackSub performs a linear solve for A x = b, given a decomposition for A. It does this efficiently by doing two simple back substitutions: one for L and one for U. The use of LUDecomp with LUBackSub to solve a system of equations is demonstrated in the following code fragment:

UT_Matrix A(1, n, 1, n);
...
UT_MatrixSolver solver;
UT_Permutation P(1, n);
fpreal determinant_sign = 0;
...
UT_Matrix LU(1, n, 1, n);
LU = A;
if( solver.LUDecomp(LU, P, determinant_sign) != 1 )
{
    // Failed to create a LU-decomposition
    return;
}
... 
UT_Vector b;
UT_Vector x;
...
// Solve A x = b
x = b;
solver.LUBackSub(LU, P, x);

Note that the LUDecomp method stores the resulting decomposition in the matrix that is passed in. For the sake of clarity, the above code fragment copies A over into a separate variable LU before calling LUDecomp. Similarly, the method LUBackSub stores the result of the linear solve in the UT_Vector parameter that's passed in.

Sparse Higher-Dimensional Matrices

The UT_SparseMatrix is used heavily in the physics simulation code. Like UT_Matrix, this class represents higher-dimensional matrices. However, UT_SparseMatrix, is optimized for the case that the vast majority of matrix elements is zero. UT_SparseMatrix only stores the nonzero elements.

The constructor for UT_SparseMatrix creates a matrix that is zero everywhere (no actual elements are stored in this case). The addToElement method can be used to initialize the matrix's elements. The following code fragment constructs a sparse matrix that can be used to perform a backward Euler step for a system of springs between n point masses in 3D:

... // Initialize array springs
UT_SparseMatrix A(3 * n, 3 * n);

for( int p = 0; p < n; ++p )
for( int c = 0; c < 3; ++c )
{
    A.addToElement((3 * p) + c, (3 * p) + c, 1);
}

for( int s = 0; s < num_springs; ++s )
{
    UT_Matrix3 dfdr[2][2]; 
    ... // Store positional derivatives of spring forces in dfdr

    for( int pi0 = 0; pi0 < 2; ++pi0 )
    for( int pi1 = 0; pi1 < 2; ++pi1 )
    for( int c0 = 0; c0 < 3; ++c0 )
    for( int c1 = 0; c1 < 3; ++c1 )
    {
        const int p0(springs[s].myPS[pi0]); 
        const int p1(springs[s].myPS[pi1]); 

        A.addToElement(
            (3 * p0) + c0, 
            (3 * p1) + c1, 
            -((h * h) / point_mass[p0]) * dfdr[pi0][pi1](c0, c1)
        );
    }
}    

Note that in this code fragment, addElement is called twice for elements on the diagonal. The result will be that the diagonal contains the sum of the values that were passed to addElement.

A UT_SparseMatrix can be multiplied with a column vector of type UT_Vector like this:

UT_SparseMatrix A(m, n);
UT_Vector x(0, n - 1);
... // Initialize A and x

UT_Vector b(0, m - 1);

// b = A x
A.multVec(x, b);

UT_MatrixSolver.h contains numerical algorithms that approximate solutions for large sparse linear systems A x = b. The UT_MatrixIterSolver class implements the preconditioned conjugate gradient method (PCG). There are various methods with the name PCG in this class. Of these, the functor-based PCG interface is the cleanest and the most flexible one. Its use is demonstrated in the following code fragment:

UT_SparseMatrix A(n, n);
UT_Vector b(0, n - 1);
... // Initialize A and b

UT_Functor2<void, const UT_Vector&, UT_Vector&> multiply_A(
    &A, &UT_SparseMatrix::multVec
);
const DiagonalPreconditioner diagonal_preconditioner(A);
UT_Functor2<void, const UT_Vector&, UT_Vector&>
    preconditioner_A(diagonal_preconditioner);

const IterationTester iteration_tester;
UT_Functor2<bool, int, const UT_Vector&> keep_iterating(iteration_tester);

// Find approximate solution x for A x = b
UT_Vector x(0, n - 1);
UT_MatrixIterSolver::PCG(
    x, b, 
    multiply_A, 
    preconditioner_A,
    keep_iterating
);
// Now A x = b (approximately)

Here multiply_A is a functor that multiplies A with a vector; it simply invokes multVec on A. The functor preconditioner_A is a preconditioner for A. The diagonal preconditioner could be implemented as follows:

class DiagonalPreconditioner
{
public:
    explicit DiagonalPreconditioner(const UT_SparseMatrix& A) :
        myID(new UT_Vector(0, A.getNumCols() - 1))
    {
        UT_Vector& id(*myID);
        A.extractDiagonal(id);
        // This asserts that A's diagonal is nonzero
        for( int i = 0; i < id.length(); ++i )
            id(i) = 1 / id(i);
    }

    void operator()(const UT_Vector& b, UT_Vector& x) const
    {
        x = b;
        x *= *myID;
    }

private:
    // inverse diagonal
    UT_SharedPtr<UT_Vector> myID;
};

The UT_Vector inside this is wrapped in a UT_SharedPtr because the DiagonalPreconditioner object is copied inside UT_Functor2. The functor keep_iterating is invoked by PCG to decide whether PCG should keep iterating, based on the current number of iterations and the residual r = A x - b. An implementation could look like this:

class IterationTester
{
public:
    IterationTester() {}

    bool operator()(int iteration, const UT_Vector& r) const
    {
        if( iteration > myMaxIterations )
            return false;
        return ( r.norm2() > myThreshold );
    }
private:
    ... // data members
};