SimTK Simmatrix™

A SimTK Toolset for efficient manipulation

of vectors and matrices in C++
Michael Sherman
Version 2.1, January, 2010
Abstract
We describe the goals and design decision behind Simmatrix, the SimTK matrix and
linear algebra library (toolset) for C++ programmers, and provide reference informa-
tion for using it. The idea is to provide the power, naturalness, and flexibility of Mat-
lab™ from within a C++ program, but with maximal performance and convenient
interoperability with numerical libraries and custom code which may already exist in
various languages, including C and FORTRAN.

1 Purpose of this document ......................................... 1 6.7 Operator reference TBD ......................................... 16

2 Goals ........................................................................ 1 Acknowledgments ........................................................ 16
2.1 Speed ........................................................................1 References .................................................................... 16
2.2 Accuracy ...................................................................2
2.3 Expressive power ......................................................2
2.4 API stability ..............................................................3
3 Design issues ............................................................ 3 1 Purpose of this document
3.1 Naming .....................................................................3
3.2 Indexing ....................................................................4
To describe the goals, design issues, theory, im-
3.3 More design issues TBD ...........................................4 plementation and use of the SimTK Simmatrix C++
4 Scalars ...................................................................... 4 toolset. Simmatrix is a part of the basic infrastruc-
4.1 Precision types ..........................................................4 ture for use of the SimTK Core tools. It is delivered
4.2 Numbers ...................................................................4 as part of the SimTKcore project on [Link]
4.3 Scalar types ...............................................................5
4.4 Scalar summary ........................................................5
([Link] and developed
5 Composite numerical types (fixed-size vectors & within the SimTKcommon project there
matrices)......................................................................... 5 ([Link]
5.1 Memory layout of CNTs; packed CNTs ...................6
5.1.1 CNT packing vs. compiler packing ............................ 7
5.2 Construction and assignment of CNTs .....................7 2 Goals
5.3 Operators on CNTs ...................................................8 Here is what we hope to achieve with Simmatrix.
5.3.1 Element access ............................................................ 9
5.3.2 Arithmetic ................................................................... 9
5.4 Summary of CNTs ..................................................11 2.1 Speed
6 Types for linear algebra.......................................... 11
6.1 Large Vector and Matrix types ...............................11 Perhaps it seems crass to start with this topic, but
6.2 Available storage types TBD ..................................12 computer simulations are dominated in practice by
6.3 Matrix characteristics ..............................................13 performance considerations. As a result, few practi-
6.3.1 Matrix character commitments ................................. 13
6.3.2 Element type ............................................................. 13
tioners will make use of a facility, however lovely,
6.3.3 Outline....................................................................... 13 that slows down their programs. A rule of thumb
6.3.4 Size ............................................................................ 14 among computational scientists is that anything
6.3.5 Structure .................................................................... 14
6.3.6 Conditioning ............................................................. 14 more than 20% overhead will begin to affect adop-
6.3.7 Sparsity ..................................................................... 14 tion by high-end users.
6.3.8 Storage formats ......................................................... 14
6.4 Matrix views ...........................................................15 When designing a system for dealing with matrices,
6.4.1 Element filters ........................................................... 15 there are two completely different and mutually
6.5 Factorizations..........................................................15
incompatible performance issues that must be ad-
6.6 Available factorizations TBD .................................16
dressed. Computations with small vectors and ma-
© 2005-10 Stanford University and Michael Sherman 1
Simbios, SimTK, Simmatrix, Simmath, and Simbody are trademarks of Stanford University.
trices, such as the 3-element ones needed for ment in existence. In fact, before you proceed fur-
representing points and orientations in 3d, are dom- ther you may wish to consider whether your needs
inated by overhead because each calculation is so might be better met by using Matlab itself. Our
small. Computations with large matrices, on the facility is not intended as a replacement for Mat-
other hand, are dominated by large-scale repetitive lab—we are interested in supporting programmers
computations and memory access patterns. Perfor- who need to work in C++ for one well-considered
mance of large computations can be dramatically reason or another. For those programmers, we
improved by a constant-time overhead spent deter- would like to provide a Matlab-like capability ac-
mining the optimal execution strategy. cessible in a natural way from within a C++ pro-
gram. We want to minimize the mental gymnastics
To deal with these conflicting concerns successful-
required to go from a mathematical statement of a
ly, we have adopted the commonly-used approach
computation to its functional implementation.
of building what are essentially two completely
Among other things, that means that real and com-
independent facilities, one suited for small objects,
plex numbers should be supported equally so that
and one for large. In the former, we do anything
the occasional appearance of a complex result does
necessary to avoid overhead. In the latter, we do
not present an insurmountable disaster.
whatever it takes to optimize computation and
memory access patterns. The facilities are then Like Matlab, we take the mathematicians’ view-
designed to work smoothly together, without hav- point that a vector is a vertical object (a column)
ing to compromise performance. A purely interpre- that is distinct from a horizontal object of the same
tive environment like Matlab can optimize only the dimension (a row or covector). We are able to en-
performance of the large computations, while con- force that distinction at compile time with no over-
tinuing to incur the same overhead on small ones. head and to do so provides significant benefits. A
matrix is then seen conceptually as a collection of
2.2 Accuracy column vectors, making ours fully compatible with
both the mathematical treatment of matrices and the
Performing linear algebra correctly on a computer
computational treatment as embodied in today’s
is a completely different field from linear algebra
best linear algebra software. However, the facility
as taught in mathematics classes. The finite preci-
is flexible enough so that anyone who prefers to
sion of computers is the dominant issue in perform-
treat matrices as collections of rows (as is common
ing correct computation, while that issue is
among computer scientists) can do so without pe-
completely irrelevant in mathematics. The best
nalty, although we would advise computer scien-
strategy here is to stand on the shoulders of the
tists not to do this reflexively—why not follow
giants who have devoted their careers to computa-
mathematical conventions for mathematical ob-
tional methods for linear algebra. In current prac-
jects?
tice, that means performing linear algebra using
software produced by the decades-long U.S. De- As mentioned above, we have a goal of zero over-
partment of Energy (DOE) effort to produce relia- head for dealing with small objects, something
ble linear algebra, as embodied in LAPACK and which is not attainable in an interpreted system.
related public domain facilities. This is the tech- C++ is one of the very few languages in which such
nique adopted by the extremely successful Matlab a facility can be attempted. Its inheritance, tem-
package, and we use it here for our large-matrix plates, operator overloading, inline functions, and
facility. We cannot use it directly for the small naughty loopholes permitting direct hardware
matrix facility when overhead concerns dominate access when necessary offer the opportunity to
everything else, but it is always available even for build extensions which provide elegant, type-safe
small systems in those cases where highest accura- abstractions with no runtime overhead whatsoever.
cy is paramount. We are able to provide zero-overhead operators for
basic matrix operations like transpose, arithmetic,
2.3 Expressive power extraction (without copying) of elements and ele-
ment subsets such as rows, columns, submatrices,
The goal here is best stated in terms of Matlab, by
diagonals, real or imaginary parts, etc. As a con-
far the most successful matrix handling environ-
crete example, matrix transpose can be seen as a

2
change in point of view rather than a physical oper- independent. Vector length is a critical attribute in
ation and C++ is powerful enough to support that programs which use the small-vector system; for
concept so that our (Hermitian) transpose ―opera- example, 3d simulation code is written using 3-
tor‖ performs no computation or memory opera- vectors and 3x3 matrices. Thus it is appropriate for
tions at all. small-vector type names to include lengths. On the
other hand, for programs written to deal with large
2.4 API stability matrices, the specific size is an unimportant run-
time specification and must not be embedded in the
An interface (or API) must satisfy very restrictive
type name.
criteria compared to general programming. Stabili-
ty is probably the primary one—an interface should To address these issues, we have adopted what we
be extremely stable once defined because many hope is a suggestive convention of using short
programs will depend on it for their correct func- names and abbreviations for short, fixed-size vec-
tioning. tors and longer ones for the large objects. Length
always appears in the short type names and never in
SimTK sets the bar higher by promising binary
compatibility. That means that software that de- the long ones. So a 3-element vector is a Vec<3>
pends on the SimTK Core interface can take advan- (abbreviated Vec3) while a long but size-
tage of new releases without being recompiled or unspecified vector is called a Vector. Similarly,
relinked in the case of dynamically linked library Mat<3,4> (abbreviated Mat34) is a 3x4 fixed-size
upgrades. This level of stability requires that ob- matrix at default precision, while Matrix is an
jects which appear in the interface must have either mxn adjustable-size matrix at default precision. The
very simple, obviously permanent implementations, full set of names is provided later in this document.
or opaque implementations which protect client One minor annoyance in C++ is that the standard
software from inevitable changes to the implemen- does not guarantee that a templatized type can ap-
tation. pear without a template argument list, even if there
All symbols introduced into a user program by are defaults for all the arguments. Thus if we want
Simmatrix are in the SimTK namespace, except for to allow Vector<Complex> some compilers will
#define macros (used sparingly) which are pre- insist on the default invocation being written Vec-
fixed by ―SimTK_‖ instead. So if we mention a tor<>, which we deem too ugly for human con-
symbol named MyType below, its actual name is sumption. We want the simple name Vector to be
SimTK::MyType. used for the by-far-most-common case, Vec-
tor<Real>. To address this, the templatized ver-
sion of the large vector and matrix types are written
3 Design issues with a trailing underscore, like Vec-
Here we discuss our thoughts on some of the choic- tor_<Complex>. Then Vector is a typedef to
es and dilemmas that are forced on anyone building Vector_<Real>, the default. This is not necessary
a system like this, and the resulting design deci- for small vector and matrix types since they always
sions we made. have at least a length argument supplied. Typical
declarations look like this:
3.1 Naming
Because there are two separate facilities, we are at
times in need of several names for similar concepts,
such as ―matrix.‖ A great deal of information must
be embedded in the types of the small matrix ob-
jects, which in C++ risks either a proliferation of
obscure type names or frequent use of templates.
Further, there is a need to permit programs to be
written as generically as is appropriate. For exam-
ple, a programmer should not need to specify preci-
sion in a type name if the algorithm is precision-

3
 Vec<3> v; // a 3-vector of Reals
Vec3 w; // same thing, using abbr. 4 Scalars
Vector b,x; // vectors of reals We start by describing the Simmatrix scalar types,
Matrix M; // mxn matrix of reals
Matrix_<Complex> C; // mxn matrix of complex that is, types which represent a single floating point
Vector_<Vec3> v3;// big vector of 3-vecs number (real or complex). We take a somewhat
novel view of these types which allows us to
// This type is a 2-element vector whose elements
// are 3-vectors. Memory layout and computational achieve some dramatic performance improvements
// efficiency are identical to Vec<6>.
when we begin to aggregate them into matrices.
typedef Vec<2,Vec3> SpatialVec; The Simmatrix scalar types comprise three levels
we call precisions, numbers, and scalars.

3.2 Indexing 4.1 Precision types

One of the thorniest issues to decide is how to treat
The precision types are the built-in C++ floating
indices, specifically, what is the index of the first
point types float, double, and long double
element? Scientific programmers used to
(called single, double, and quad or real*4,
FORTRAN, Matlab, and general mathematical con-
real*8, and real*16 in FORTRAN). These con-
ventions expect the first element to have index ―1.‖
C programmers expect indexing to begin with ―0.‖ vey the level of accuracy required in computations
Many packages address this by making the index- and stored values. On all currently supported plat-
ing offset a user choice. We have done this in past forms, float is 4 bytes, and double is 8 bytes.
designs and found it unsatisfying and extremely Depending on the compiler being used, long
prone to induce errors where some programmers double may be the same as double, or it may be
(typically at the upper levels of the code) use 1- 10 or 16 bytes. All SimTK-supported computers
based indexing while others use 0-based. The re- must adhere to the IEEE 754 standard for floating
sulting awkwardness and uncertainty produces point arithmetic, and SimTK-compliant code may
either subtle off-by-one errors, or unnecessary co- depend on that fact, for example by depending on
pying as responsible programmers move data into the existence of NaN (not-a-number) and Infinity.
compatibly-indexed vessels to avoid errors. Under compile-time control, one of the precisions
So we believe that a single indexing scheme must is chosen as the default precision and given the type
be chosen and used exclusively. Following the name Real (reminder: this is really
reasoning of the VXL design team (see SimTK::Real). We expect most SimTK code to
[Link] be written in terms of Real rather than the explicit
e/book_6.html) we have opted for 0-based indexing precisions, so that it may be compiled at different
as ―least weird‖ in a C++ numerical package. This precisions without changing the code. Our default
does not preclude a FORTRAN-compatible 1-based for this default is that Real is equivalent to
interface to the facility, but it prevents any leakage double, that is, an 8-byte floating point value with
of incompatible indices into the guts of the pack- approximately 16 decimal digits of precision and an
age. This will not please everyone, but in past ef- exponent range of about ±300 (base 10).
forts to do so many packages have ended up
Contributors to SimTK code should avoid writing
pleasing no one.
precision-specific code whenever possible. In addi-
tion to our predefined default-precision types, C++
3.3 More design issues TBD provides the standard template type
Topics: choice of indexing operators, subvector and std::numeric_limits, and we provide all the
submatrix operations, layout of symmetric matric- additional tools needed occasionally for writing
es, mixed-precision operations, choice of transpose precision-independent code, such as generic-
operator, use of non-conforming operands, why precision constants and convergence limits.
can’t large matrix and vector types be used as ele-
ment types, packing of data in memory, treatment 4.2 Numbers
of scalar assignment.
From the above precision types, we construct our
numbers, each of which can exist in all three sup-

4
ported precisions. The number types are the stan- 4.4 Scalar summary
dard real and complex numbers, and a novel ―fla-
The above defines a set of exactly 18 scalar types:
vor‖ of complex number called a conjugate, which
three kinds of numbers in each of three precisions,
is not typically used in user programs but is impor-
in normal or negated form. Despite the use of tem-
tant for the efficient implementation of matrix op-
plates, this is not a user-extendable set.
erations.*
Stated as a grammar, scalars are defined like this:
Each of the three precision types is also a real
number type, with the SimTK type Real prede- scalar ::= number | negator<number>
fined equivalent to one of those as described above. number ::= standard | conjugate
The C++ standard template type std::complex, standard ::= real | complex
conjugate ::= Conjugate | conjugate<precision>
specialized as complex<float>, com-
complex ::= Complex | complex<precision>
plex<double>, and complex<long double>
real ::= Real | precision
serve as our complex numbers, with SimTK de- precision ::= float | double | long double
fault-precision type Complex predefined as com-
plex<Real>.† The three conjugate types (from the That completes our discussion of the scalar types.
SimTK namespace) are conjugate<float>, Next we’ll look at how these can be used to con-
conjugate<double>, and conju- struct the much more interesting composite numeri-
gate<long double>, with default-precision type
cal types.
Conjugate predefined as
conjugate<Real>. 5 Composite numerical
types (fixed-size vectors &
4.3 Scalar types matrices)
Only Real and Complex, and vectors and matrices
SimTK defines the following composite numerical
defined in terms of them, will appear in typical user
types (CNTs), built up recursively from scalars and
programs. However, our complete set of scalar
other composite numerical types. As a consequence
numeric types consists of the nine number types
of our definition of ―scalar‖ above, all composite
described above, plus a templatized adaptor class
numerical types support both real and complex
negator<number>, which may be applied to any
arithmetic equally.
number type to create a new type whose memory
representation is unchanged but whose value is to CNTs can be scalars as described above, or small
be interpreted with the opposite sign. Like conju- fixed-size vectors and matrices whose elements are
gate, negator is not expected to appear in user scalars or other CNTs. For example, one may easily
programs but permits efficient implementation of define a 2-element vector whose elements are ordi-
some matrix operations, in particular Hermitian nary 3-element vectors.‡ CNTs are characterized by
transpose of a complex matrix. The imaginary part zero overhead and minimal storage requirements.
of a conjugate number has type For example, an mxn matrix of scalars is stored as
negator<real> for the appropriate precision real. exactly mn consecutive scalars (except when oth-
erwise requested), with no flags, counters, heap
*
pointers, or other sources of inefficiency. The exact
Conjugate types have the same representation as com- storage layout in memory is part of the definition of
plex types, but the imaginary part is interpreted with
these types (and is machine-, operating system-,
opposite sign.
and compiler-independent), so they may be recast
†
The one difference between our complex type and the to native machine types or other CNTs with pre-
C++ standard is that we do not initialize unused values dictable results. This implies, for example, that a
in a Release build. That is, we treat complex identically
to real in this regard, while the C++ standard specifies
‡
that complex variables are initialized to (0,0). We feel That is very useful as the representation of a Spatial
that is an inappropriate default for large, complex ma- Vector which collapses both rotational and translational
trices where avoiding unnecessary memory accesses is effects into a single quantity in the Spatial Operator
of primary concern. Algebra used by SimTK’s Simbody™ package.

5
symmetric matrix has a different type from a gen- SymMat<m> A small, fixed-size mxm symmetric
eral one, so that the appropriate access pattern can SymMat<m,C> (Hermitian if complex) matrix of
SymMat<m,C,rs> composite numerical type C (default
be determined at compile time. CNTs make exten- Real) with optional element-to-
sive use of C++ templates and inline functions, and element spacing (default 1). Only the
the fact that type casting is free, to permit compile- SymMat22 elements of the diagonal and lower
time optimization into the optimal set of machine SymMat33 triangle are stored.
instructions. SymMat44 These typedefs are provided, with
… SymMat99 SymMat33 ≡ SymMat<3>, etc.
The table below presents the available Composite
Numerical Types, or more precisely the templates
available for constructing CNTs. Note that prede-
Note that the short Vec and Mat types are them-
fined abbreviations are available for certain com-
selves Composite Numerical Types and can thus be
mon combinations of template arguments (up to 9
composed recursively. That is, it is reasonable to
rows and columns). We provide those because we
know from experience that programmers will de- have a 2x2 matrix of 3x3 matrices and in fact this
can be quite useful. This can be declared*
fine them to reduce the ugly ―<>‖ clutter created by Mat<2,2,Mat<3,3> >
the C++ template syntax, and we would like to or more pleasantly using a predefined typedef
encourage a consistent set of common abbrevia- Mat<2,2,Mat33>
tions. These abbreviations are not distinct types Note that these types are exactly equivalent and
(they are typedefs), so may be freely intermingled thus interchangeable, and that the resulting type is
with the spelled-out types. For example, a Vec3 itself a CNT.
can be passed to a routine written to take a Vec<3>
argument. 5.1 Memory layout of CNTs;
Type Description packed CNTs
Real A floating point number at default The default layout for any CNT is to pack the ele-
Complex precision, either real or complex. ments into the least amount of consecutive storage
Vec<m> A short, fixed-length column vector of as logically required to hold the CNT’s value. We
Vec<m,C> m elements of composite numerical refer to a CNT stored this way as a packed CNT. In
Vec<m,C,stride> type C (default Real) with optional addition, each of the CNT templates provides ar-
element-to-element stride (default 1).
guments which can be used to specify regular gaps
Typedef abbreviations are provided as
Vec2 Vec3 Vec4 in the storage layout, where the gap size is always
shown, with Vec3 ≡ Vec<3>, etc.
Vec5 … Vec9 an integer multiple of the storage requirement of
Row<n> A short, fixed-length row vector of n the CNT’s element type. For example, the Vec and
Row<n,C> elements of composite numerical type Row templates allow specification of a stride, which
Row<n,C,stride> C (default Real) with optional ele- gives the spacing between consecutive elements in
ment-to-element stride (default 1).
Rows are not typically used in user terms of those elements. So a packed Vec or Row
Row2 Row3 Row4 programs, but occur as intermediate has stride 1, which is the default.
Row5 … Row9 results in expressions.
Non-packed CNTs exist to facilitate reinterpreta-
Mat<m,n> A small, fixed-size matrix of m rows
Mat<m,n,C> and n columns of composite numerical tion of existing data in terms of CNTs. For exam-
Mat<m,n,C,cs,rs> type C (default Real) with optional ple, if one has a Mat<4,3> stored as three
column-to-column spacing cs (default consecutive packed Vec<4> columns, the rows can
m), and row-to-row spacing rs (default be viewed as a 3-element Row CNT with a stride of
Mat22 Mat33 … 1).
4 (that is, four Real elements), which could be
Mat66 … Mat76 These typedefs are provided, with
… Mat89 Mat99 specified as Row<3,Real,4>. However, such dec-
Mat33 ≡ Mat<3,3>, etc.
larations do not normally appear in user programs;

*
C++ unfortunately requires the extra space for nested
templates, to avoid confusion with the operator ―>>‖.

6
instead, they are the hidden return types of methods same elements. So while you can always cast a
and operators which select out portions of an exist- CNT into a C++ scalar array with predictable re-
ing object. In this case, the row index operator m[i] sults (for any of the supported scalar types includ-
acting on a Mat<4,3> matrix m returns the ith row ing complex), you cannot safely cast to a C++ array
of m as a 3-element Row with stride 4, meaning it of composite CNT types; cast to a CNT Vec of
references the elements of m without copying, and those types instead.
can serve as an lvalue (target of an assignment)
which alters the appropriate elements of m. Because 5.2 Construction and assignment
this type has the same semantics as any other of CNTs
Row<3>, most users will never need to think about All CNTs define a default constructor. In Debug
how it is implemented. mode (that is, when the C++ standard NDEBUG pre-
When we discuss the large-matrix facility below, processor symbol is not defined), the default con-
the distinction between packed and non-packed structor initializes all elements to NaN. In Release
CNTs is somewhat more significant, since only mode (i.e., NDEBUG is defined), all elements are left
packed CNTs can serve as element types for the uninitialized, so that declaring CNTs, or arrays of
large Vector and Matrix classes. CNTs, has no cost if the variables are not used.
Constructors are also available for initializing data
5.1.1 CNT packing vs. compiler pack- elements from individual element values, or by
ing copying compatible CNTs. Initialization values can
Some C++ compilers attempt to improve execution be provided in the constructor or via a pointer (or C
speed by ―rounding up‖ memory requirements to a array) to values of the appropriate type.
size which is particularly efficient for the targeted
hardware. For example, a class containing three 4- Assignment operators are available for copying one
CNT to another, and for setting single elements,
byte float values (e.g., Vec<3,float>) might be
subvectors and submatrices.
allocated 16 bytes rather than 12, by some compi-
lers on some machines, with some compile-time One important convention we follow, which is
options. That means that C++ arrays of such ob- different than that of most similar systems, is the
jects would contain 4-byte gaps between the ele- treatment of scalar assignment. We follow this
ments. convention: (1) when a scalar s is assigned to a
vector, every element of the vector is set to s (this
Because a great deal of our performance comes
is the typical convention), and (2) when a scalar s is
from the ability to change our point of view (i.e.,
assigned to a matrix, the diagonal elements of the
recast) rather than copy or compute, we depend on
matrix are set to s while the off-diagonals are set to
predictable storage layouts for our classes. To en-
zero. Examples:
sure that, we define packed CNTs in terms of the
storage requirements of their underlying scalar Vec3 v;
Mat22 m;
types, which are packed the same way by all com- Vector b(10); // initial size 10 reals
pilers on all machines. We guarantee, for example, Matrix M(20,10); // initial size 20x10
that a CNT Vec<2,Vec3> can be recast to a v=0; // v=(0,0,0)
Vec<6> with the obvious interpretation and memo- v=3; // v=(3,3,3)
m=0; // m=( 0,0 )
ry layout identical to a C++ array Real[6]. This // ( 0,0 )
would not necessarily be possible if the CNT stored m=1; // m=( 1,0 )
// ( 0,1 )
the two Vec3’s in a C++ array Vec3[2] since the b=0; // b=10 zeroes
compiler might choose to allocate space for two M=1; // M=0, except M(i,i)=1, 0<=i<10
Vec4’s instead!
This convention is especially apt for matrices, be-
The key point to remember is that CNTs are packed cause the matrix resulting from such a scalar as-
internally as arrays of scalars, so that a composite signment ―acts like‖ that scalar. That is, if you
CNT whose elements are also composite CNTs multiply by this matrix the result is identical to a
may occupy less storage than a C++ array of those scalar multiply by the original scalar. Two impor-

7
tant special cases are: (1) setting a matrix to the
scalar ―1‖ results in the multiplicative identity ma-
trix of that shape, and (2) setting a matrix to the
scalar ―0‖ results in the additive identity matrix of
that shape. In general, in any operation involving a
scalar s and a Matrix or Mat, the scalar is treated
as if it were a conforming matrix whose main di-
agonal consists of all s’s with all other elements
zero. So Matrix m += s will result in s being add-
ed to m’s diagonal, which is what would happen if s
were replaced by diag(s) of the same dimension as
m. m-=1 thus subtracts an identity matrix from m,
without touching any of the off-diagonal elements.
Note that for multiply and divide this convention
yields the ordinary scalar multiply and divide oper-
ations: e.g., m*s (=m*diag(s)) multiplies every
element of m by s.

5.3 Operators on CNTs

There are numerous operators available which act
on CNTs, for element access, computation, reinter-
pretation of the data, and obtaining information
about the type and its contents.

8
5.3.1 Element access
These operators provide access to individual elements, or subsets of elements, of composite numerical types,
where we use letters to indicate types: s=scalar, e=element (of whatever CNT), v=Vec, r=Row, m=Mat,
sy=SymMat. i,j are integer indices, with i a row index and j a column index. An ―Lvalue‖ can appear on the
left hand side of an assignment statement, with the result affecting the original element values.

Operator Applied Meaning Lvalue? Cost Notes

to
v[i] v(i) Vec select ith (jth) element yes native array index all indexing is 0-based
r[j] r(j) Row
m[i][j] m(i,j) Mat obtain i,j element of m. yes same as native
SymMat Only diag & lower matrix index for
sy[i][j] sy(i,j) triangle of SymMat; Mat; extra integer
i.e., i ≥ j. operations for
SymMat.
m[i] [Link](i) Mat obtain ith row or jth yes native array index Size and spacing are taken
m(j) [Link](j) column of m as Row or from Mat; typically col-
Vec, resp. umns are packed while
Rows have stride>1.
[Link]() Mat obtain diagonal as a yes zero For rectangular Mat<m,n>,
[Link]() SymMat Vec result has dimension
min(m,n).
sy[i] [Link](i) SymMat obtain ith row or jth no must copy ele- Return type is packed
sy(j) [Link](j) column of SymMat as a ments to temporary regardless of original
Row or Vec Row or Vec; avoid SymMat spacing.
if possible
[Link]<m>(i) Vec return a reference to a get – no native array index Element type and stride are
[Link]<m>(i) Vec<m> whose 0th upd – yes the same as the original
element is v’s ith ele- vector.
ment
[Link]<n>(j) Row return a reference to a get – no native array index Element type and stride are
[Link]<n>(j) Row<n> whose 0th upd – yes the same as the original
element is r’s jth ele- row.
ment
[Link]<m,n>(i,j) Mat return a reference to a get – no native array index Element type and spacing
[Link]<m,n>(i,j) Mat<m,n> whose 0,0 upd – yes are the same as the original
element is m’s i,j ele- matrix.
ment
[Link]<m>(i,j) Mat return reference to a get – no native array index Element type and stride are
[Link]<m>(i,j) SymMat Vec<m> or Row<n> upd – yes the same as the original
whose 0th element is vector.
[Link]<n>(i,j) m’s i,j element. Only
[Link]<n>(i,j) strictly lower triangle
(upd also available) can be referenced for
SymMat; that is, i > j.

5.3.2 Arithmetic
The expected arithmetic operators are overloaded for use with CNTs, plus probably some unexpected ones.
This includes add, subtract, matrix multiply, and divide for conforming objects and scalars, cross product for

9
2- and 3- vectors, and the usual C arithmetic assignment operators like ―+=‖. Behavior for CNTs with scalar
elements are as expected; behavior for CNTs with composite elements are defined analogously and generally
work well, but most users will not have a well-developed intuition for those objects at first.

Operator Applied Meaning Lvalue? Cost Notes

to
~ (transpose) Any CNT Transpose (Hermitian yes zero (implemented Has no effect on symmetric
transpose if elements as a cast) matrices.
are complex). Acts as ―conjugate‖ operator
on complex scalars.

+-*/ Any CNT matrix arithmetic no same as explicit conformant matrices always
+= -= *= /= code work; some non-conformant
Dot and outer Row*Vec is a dot operations are useful also
s=dot(v,w) work for any product (scalar result); Global methods
combination Vec*Row is an outer have the same dot() uses the Hermitian
m=outer(v,w) transpose of v’s elements times
of Vec & product (Mat result). performance as the w’s unchanged elements,
Row if corresponding regardless of what is a Row or
lengths match. Global methods dot() operators. Vec.
and outer() are
provided as an explicit Similarly, outer()uses v
alternative. unchanged and the Hermitian
transpose of w.

% (cross product) 2- and 3- Returns cross product. no Same as explicit If both arguments are Vec3,
element Vec Result is a scalar for 2- code (3 flops for 2- result is Vec3. If either is
z=cross(v,w) and Row. element cross product, element, 9 flops Row3, then result is Row3.
and a 3-vector for 3- for 3-element). 2-element cross product can be
element cross product. understood as 3-element where
Global method and a zero z component has been
Global method operator provide added and the final result is the
cross() is provided the same perfor- z component of the result.
as an alternative. mance.

m=crossMat(v) 2- and 3- Returns the matrix no same as explicit The result is a 3x3 skew-
element Vec which acts as a cross code (a few nega- symmetric matrix in the 3d
and Row. product operator. That tions and copies) case, a 2-element Row in the
2d case. The same result is
is, mw==vxw.
produced regardless of whether
the argument is a Vec or Row.

10
5.4 Summary of CNTs
Type Description
Stated loosely as a grammar, we build CNTs recur-
sively like this: Vector An arbitrary-length column of Real
Vector_<C> values, or of values of packed compo-
CNT ::= scalar site numerical type C (e.g.,
| composite<size [,CNT [,packing] ]> Vector_<Complex> or
composite ::= Vec | Row | Mat | SymMat Vector_<Mat<2,2,Mat33> >).
size ::= nrow [,ncol] RowVector Same as Vector but horizontal.
packing ::= stride | colSpacing,rowSpacing RowVector_<C> Usually not used explicitly in code, but
is the type of a Matrix row or
The unbolded terminals nrow, ncol, stride,
Vector transpose.
colSpacing, and rowSpacing are integers, or
compile-time expressions that evaluate to integers. Matrix An arbitrary-size, two dimensional
Matrix_<C> matrix of Real values, or of values of
This grammar permits unlimited nesting of these packed composite numerical type C.
constructs, and the implementation does work that
way, but we would counsel restraint here and note
that there is unlikely to be much utility (or clarity) Standard linear algebra operations, matrix decom-
beyond two or three levels deep. positions, and interconversions with composite
numerical types are provided. Note that SimTK
6 Types for linear algebra Vector and Matrix are not themselves composite
numerical types and may not be composed recur-
[This part of the document is very sparse at the sively. A Vector may be considered an mx1 Ma-
moment (no pun intended).]
trix, and a RowVector a 1xn Matrix when
convenient. Thus the discussion below which refers
6.1 Large Vector and Matrix types to matrices applies to Vector and RowVector as
The ―zero overhead‖ requirement on the Composite well.
Numerical Types limits their flexibility. For larger
Unlike the Composite Numerical Types, very little
vectors and matrices, some constant-time overhead
is encoded in the type here—only the basic shape
is acceptable since we expect time to be dominated
by floating point calculations and memory accesses outline (1 or 2d object) and the element type. Di-
done on the (large) operands. In fact, this overhead mensions, spacing, and internal data layout are
is desirable since it is used to set up optimal large- determined at runtime. This provides a great deal of
scale operations which can then be performed at useful flexibility but imposes a constant-time cost
for every operation.
machine speeds. The basic types, and general beha-
vior, are modeled after the very successful Matlab SimTK provides 0-based indexing using the []
system with the expectation (but not requirement) operator. If the Matrix is modifiable (non-const)
of LAPACK and BLAS style implementation. We then the indexed element can be modified and that
assume that these objects will be very large and the change affects the contents of the object. The []
classes are carefully designed to avoid unnecessary operator applied to a Matrix returns a row, which
data copying and memory references. may in turn be indexed to obtain an element in C
The underlying element type stored in our large style. SimTK also permits indexing using round
matrix objects can be any scalar type or other brackets () yielding identical results to [] for
packed composite numerical type (see section 5.1). Vector but selecting a column rather than a row
These elements will be packed adjacent in memory when applied to a Matrix. A two-argument round
in our large matrix objects regardless of whether bracket operator accesses a Matrix element, and
the C++ compiler would pack them that tightly unlike for CNTs, it is more efficient to use the two-
when creating its own arrays. Other data layouts are argument form here since the overhead cost is paid
available if explicitly requested, but packing of only once.
elements is always done by packing the underlying Matrix m; Vector v; …
scalars as discussed for CNTs in section 5.1.1. v[i] // ref to ith element of v, 0-based

11
 v(i) // same posite types, and that a negated scalar is still a sca-
m[i][j] // ref to i,jth element of m, 0-based lar. See section 4 for more information about
m(i,j) // same, but faster
SimTK scalars.
m[i] // ref to ith row of m, 0-based
m(j) // ref to jth column of m, 0-based Many operations with composite CNT elements
There are also operators for selecting subvectors can be performed at full speed also, provided that
and submatrices. Like the indexing operators, these there is an equivalent scalar operation. For exam-
return references into the original object, not cop- ple, scalar multiplication on composite elements is
ies. Submatrices are thus ―lvalues‖ (in C terminol- equivalent to the same operation done on the ele-
ogy) meaning that they can appear on the left hand ments’ underlying scalars.
side of an assignment. The underlying data representations for these ob-
Matrix m; Vector v; … jects are documented and stable, and in many cases
v(i,m) // ref to m-element subvector whose 0th map directly onto standard dense storage. In those
// element is v’s ith element cases a pointer to the raw data can be obtained if
m(i,j,m,n) // ref to mxn submatrix whose (0,0) necessary for performance or compatibility purpos-
// element is m’s (i,j) element es.
References of this type are called views since they
provide alternate views of the same data. They 6.2 Available storage types TBD
retain all properties of the original object except Most Simmatrix physical data layouts are designed
that they cannot be resized. They are in fact to be directly compatible with one of the LAPACK-
represented identically to the original objects in the defined layouts. The default layout assumes that an
sense that they can be used wherever a Vector or mxn matrix is stored by columns using mxn ele-
Matrix reference is expected, without memory ments. An option is to provide a ―leading dimen-
allocation or data copying. sion‖ (≥ m) so that there are regular memory gaps
The implementations of these types are opaque to a between the columns. (This usually occurs auto-
C++ program using them. That is, the header files matically when selecting a block view of part of a
define these as ―handle‖ classes which contain only larger matrix.)
a pointer to an undefined type (essentially a Our default for symmetric and triangular matrices
void*). The object referenced is an instance of a is what LAPACK calls ―conventional‖ storage; that
hidden implementation class.* This permits use of is, space is allocated for the whole matrix, but only
these classes in SimTK interfaces while preserving half the space is used. The space-saving ―packed‖
binary compatibility. It also allows use of these format is also available but can be expected to run
objects from other languages, since the void* can significantly slower for many operations. Note that
serve as a ―lowest common denominator‖ represen- many LAPACK operations pack two symmetric or
tation. triangular results into a conventional matrix; that
As in Matlab, there is a substantial performance provides both optimal speed and space.
penalty to work with Vector and Matrix objects Banded matrices are also available, always in
element-by-element; ―bulk‖ operators should al- packed storage. These can be full, symmetric, or
ways be used in performance critical code. For triangular.
some operations, matrices with structured CNT
We provide direct support for permutation matric-
elements will perform poorly compared to matrices
es, which are permuted identity matrices resulting
of scalar elements. When used properly, Simmatrix
from pivot operations required for numerical stabil-
objects containing scalar elements are capable of
ity during factoring. The underlying storage for
performing large-scale operations at full machine
these is typically just a sequence of integers defin-
speed; that is, as fast as LAPACK. Note that com-
ing how the columns or rows were pivoted.
plex and conjugate types are still scalars, not com-
(TODO: I think LAPACK has two different layouts
of pivot matrices; we’ll have to support them both.
*
This is a standard C++ design pattern usually called It should be possible to extract the underlying in-
―PIMPL‖ for ―private implementation.‖ teger array and pass it straight through to LAPACK;

12
by hiding it under a Matrix handle we don’t have 6.3.1 Matrix character commitments
to worry about whether the integer indices start at In general, a Matrix handle will be committed
zero.) with respect to some of the above attributes, and
TODO: ideas for more: scalar matrices for zero and uncommitted to the rest. A character commitment
identity and scalar*identity; sparse matrices in specifies ―minimum acceptable‖ properties for the
some DOE-compatible format(?) actual matrix referenced by the handle. For exam-
ple, a matrix handle committed to symmetric struc-
Where possible, similar storage options are availa- ture cannot be assigned to a nonsymmetric result,
ble for any element type, however factoring, pivot- but would accept a diagonal result, since every
ing and so on are only defined for scalar elements. diagonal matrix is also symmetric.

6.3 Matrix characteristics Every Matrix (and Vector and RowVector) is

committed to a particular element type, since an
A declaration like ―Matrix m;‖ declares m as an element type is required as a template argument in
uncommitted matrix handle. That is, although m has the Matrix type itself (recall that Matrix itself is
the semantics of a 0x0 matrix of reals, it has not yet an abbreviation for Matrix_<Real>). In addition,
been committed to using a particular data layout. It Vector and RowVector handles are committed to
can be used to hold the results of any matrix opera-
a particular outline: column and row, respectively.
tion, and will take on the characteristics of that
result. A subsequent use of the handle might leave Any other desired commitments must be added
it with completely different characteristics; only the explicitly before the handle is used to hold any
element type can never change. If an uncommitted data. Many of the characteristics represent catego-
handle is asked to allocate space for some data (for ries of acceptable attributes, rather than specific
example, via ―[Link](10,20);‖) it will use a ones. For example, a commitment to a square out-
dense, column oriented allocation identical to line still permits flexibility with regard to the size,
LAPACK’s conventional matrix storage format. as long as both dimensions are the same. The ma-
trix characteristics are not completely indepen-
However, handles can optionally be restricted to dent—some imply others. For example, a
narrower ranges of behavior, via commitments symmetric structure implies a square outline.
which will be discussed below. First we’ll discuss
the kinds of characteristics that a matrix can pos- Next we’ll look at the individual matrix characteris-
sess. tics one by one.
A Simmatrix Matrix, Vector or RowVector 6.3.2 Element type
object is characterized by the following seven
attributes: As mentioned above, every matrix handle is com-
mitted to a particular element type by its own tem-
1. Element type (a CNT) platized type. Only packed CNTs are permitted as
2. Outline element types.
3. Size
4. Structure Most operations require exactly matching element
5. Conditioning types among the operands, with the exception that
6. Sparsity operands which differ only in the negation or con-
7. Storage format jugation status of their underlying scalars can be
intermingled.
Collectively, we refer to a set of particular values
of these attributes as a matrix character. There are 6.3.3 Outline
two matrix characters associated with every matrix There are five possible outline attributes:
handle: the handle’s character commitment, and its
current character. The current character always 1. Rectangular (mxn)
satisfies the handle’s commitment. 2. Square (nxn)
3. Column (mx1)
4. Row (1xn)

13
 5. Scalar (1x1) 6.3.6 Conditioning
A matrix handle with no outline commitment can Matrix condition is a statement about the numerical
hold a general (rectangular) matrix, which of properties of a Matrix. It can be set as a result of an
course includes all the other outlines as well. operation, or by a knowledgeable user. Simmatrix
is entitled to rely on the correctness of these asser-
Vector handles are always committed to column tions, although it will try to check them when it is
outline, RowVector handles to row outline. possible to do so without sacrificing performance.

6.3.4 Size 1. Unknown (the default)

2. Singular (expect terrible conditioning)
Size can be variable or fixed, in one or both dimen-
3. Full rank (for rectangular, this means rank
sions. When a dimension is fixed, the matrix handle
is the same as the shortest dimension)
must commit to a particular size for that dimension.
4. Well conditioned (implies full rank)
Vector handles are always committed to exactly 5. Positive definite (symmetric only)
one column, RowVector to exactly one row. Even 6. Orthogonal
a zero-length Vector has one column, that is, it is Most operations will result in Unknown condition-
0x1, a zero-length RowVector is 1x0. ing, which will not satisfy a more restrictive com-
mitment. In the case that a matrix handle commits
6.3.5 Structure to particular conditioning, only operations which
Structure refers to an inherent mathematical (or at preserve that conditioning will be successful.
least algorithmic) property of the matrix rather than
a storage strategy. Symmetry is the clearest exam- 6.3.7 Sparsity
ple of this; it is far more significant than just a way This is a statement about the number and/or place-
to save storage and reduce operation count. ment of non-zero entries in a matrix. When we
1. General (the default) know that few entries will be non-zero, exploiting
2. Vector (can be a column vector or a row this fact can yield significant speed gains. It is par-
vector) ticularly easy and efficacious to work with a matrix
3. Symmetric (includes Hermitian); implies whose non-zero elements cluster narrowly about
Square outline; Hermitian implies real di- the main diagonal.
agonal. 1. Full (default)
4. Triangular (includes trapezoidal) 2. Banded: diagonal plus specified upper and
a. Upper/lower lower bandwidths
b. Hessenberg (triangular except for a. If symmetric or triangular struc-
one sub- or super-diagonal) ture, only one bandwidth
c. Quasi-triangular (triangular except b. Bandwidth may be left uncommit-
for 2x2 blocks on the diagonal) ted or commit to particular band-
5. Diagonal (also symmetric & triangular); width
rectangular outline is OK 3. TODO: Sparse.
6. Scalar (diagonal, and all diagonals have the
same value); doesn’t imply Scalar outline 6.3.8 Storage formats
7. Permutation (of rows or columns); implies
These refer to the physical layout of data in the
Square outline
computer’s memory. Whenever possible we at-
A ―vector‖ structure means that the mx1 or 1xn tempt to store data in a format that enables use of
matrix is mathematically a vector or covector, special high performance methods, such as those
which is subtly different from a mx1 or 1xn matrix. available in the SimTK LAPACK/BLAS implemen-
In particular certain norms are defined differently tation.
for vectors than for similarly-shaped matrices.
1. Full (default for full, symmetric & triangu-
lar matrices)

14
 a. Specify leading dimension (default nested. So a Matrix always contains at most one
is number of rows) view.
b. Symmetric, triangular, diagonal
Many views can be manipulated as efficiently as
can exist within full storage
the original data, particularly when the original is a
c. Specify upper/lower for symmetric
full-storage matrix.
& triangular
d. Specify whether unit diagonal is
6.4.1 Element filters
assumed
2. Packed (default for banded matrices, space Whenever possible, we construct a matrix view
saving for symmetric & triangular but simply by finding another high-performance de-
usually considerably slower than using full scription of the desired subset of the data. For ex-
storage) ample, a view which is a row of an ordinary dense,
a. Specify whether unit diagonal is column-oriented matrix can be represented as a
assumed ―strided‖ one-dimensional object, which is one of
the formats which can be manipulated efficiently
3. Householder product (LAPACK representa- by the BLAS routines.
tion for orthogonal matrices)
It is possible that a desired view cannot be ex-
4. Pivot array (set of integers used to pressed in one of the available high-performance
represent Permutation matrices) data descriptors. In that case the matrix supple-
5. TODO: look into use of rectangular full- ments the data descriptor with an element filter. An
packed storage (Gustavson & Wasniewski) element filter presents a logical matrix whose ele-
instead of LAPACK packed – claim is up to ments consist of a subset and/or reordering of the
20X faster. in-memory elements, done in a way that does not
map to a supported high speed format. For exam-
6. TODO: should we include a stride as a data ple, one could construct a view which picked out
storage option or is that just a view? particular elements, with arbitrary spacing and
ordering compared to the originals. These could
6.4 Matrix views appear as a contiguous Vector, for example, al-
In addition to being an owner of element data, a though the individual elements might be widely
matrix handle can also serve as a view into some- scattered. Operations on such an object are unlikely
one else’s element data. A view provides a logical to be very efficient, but in many cases the clarity of
matrix whose elements consist of a subset and/or code will matter more.
rearrangement of the elements described by the
data descriptor, sometimes augmented with a few 6.5 Factorizations
generic read-only data elements like 0, 1, and NaN. For equation solving, one may always calculate a
Views are commonly used to select blocks or di- matrix inverse and then multiply by it; however,
agonals from within a large matrix or to provide a this is the mathematician’s approach rather than the
matrix which appears to be transposed relative to computational scientist’s and at times will not yield
the original data. A Matrix view is an lvalue and acceptable results in finite precision arithmetic.
assignment to the view results in changes to the Simmatrix does allow that approach but it is not
actual elements of the original data. recommended. As a better option for one-time use,
The most commonly encountered views are those the divide operator is overloaded to allow casual
created by operators such as transpose or row, col- solution to Mx=b by writing x=b/M, which means
umn and submatrix selection. x=M-1b (or x=M+b if a pseudoinverse is neces-
sary). M can contain information about its condi-
A Matrix view is still a Matrix, and can be passed tioning and structure which permits the operator to
to any Matrix argument. Views may be made of make a reasonable choice of solution method; oth-
views, with the logically correct result, however the erwise, Simmatrix will make a conservative choice
views are combined into a single view rather than yielding good numerical results but perhaps subop-
timal performance. In any case the divide operator

15
does not actually form the inverse, but works di- Acknowledgments
rectly with the factorization which is numerically
preferable. This work was funded by the National Institutes of
Health through the NIH Roadmap for Medical Re-
For more control, or for repeated use of the same search,1 Grant U54 GM072970.
matrix while factoring only once, one must con-
struct explicit factorizations and then solve equa-
tions using the factorization directly rather than
References
using it to invert the original matrix.
1
Matrix factorizations are objects which can be used Information on the National Centers for Biomedical
similarly to matrix inverses, but with optimal nu- Computing can be obtained from
[Link]
merical accuracy. See the Simmath documentation
for information on Factor classes; here’s an exam-
ple:
Matrix M; Vector b1,b2,x1,x2; …
FactorLU f(M); // LU factorization of M
x1 = f*b1; // instead of x1=b1/M
x2 = f*b2;

These can be made to yield the best possible results

with the highest efficiency, and the Factor classes
can provide many useful methods such as rank
determination. Typically the Factor constructor
will obtain the layout and known properties from M
and then call the appropriate LAPACK routines to
perform the factorization. Options exist to allow the
Factor class to steal the original memory from the
matrix being factored.

6.6 Available factorizations TBD

Square, well conditioned matrix: LU with pivoting
Symmetric, general matrix: LLT
Symmetric, positive definite: Cholesky (LDLT?)
Rectangular: QR with pivoting, LQ
Rectangular, ill conditioned: QTZ, SVD
Symmetric and nonsymmetric eigenvalues routines
and Schur factorization
Access to the underlying factors (without copy-
ing!).
Condition number, rank determination/setting, equ-
ation solve, inverse and pseudoinverse.
How error conditions are handled.

6.7 Operator reference TBD

Basic Matlab and BLAS equivalents.

16