Overview of Matrix and Vector Operations

Contents:: Basic Linear Algebra; Advanced Functions; Submatrices, Subvectors; Speed Improvements

Definitions:

`A, B, C`	are matrices
`u, v, w`	are vectors
`i, j, k`	are integer values
`t, t1, t2`	are scalar values
`r, r1, r2`	are ranges, e.g. `range(0, 3)`
`s, s1, s2`	are slices, e.g. `slice(0, 1, 3)`

Note: A range r = range(start, count) contains all indices i with start <= i < start+count. A slice is something more general. The slice s = slice(start, stride, count) contains the indices start, start+stride, ..., start+(count-1)*stride. The stride can be any integer including 0 and -1!

Basic Linear Algebra

standard operations: addition, subtraction, multiplication by a scalar


C = A + B; C = A - B; C = -A;
w = u + v; w = u - v; w = -u;
C = t * A; C = A * t; C = A / t;
w = t * u; w = u * t; w = u / t;

computed assignements


C += A; C -= A; 
w += u; w -= u; 
C *= t; C /= t; 
w *= t; w /= t;

inner, outer and other products


t = inner_prod(u, v);
C = outer_prod(u, v);
w = prod(A, u); w = prod(u, A); w = prec_prod(A, u); w = prec_prod(u, A);
C = prod(A, B); C = prec_prod(A, B);
w = element_prod(u, v); w = element_div(u, v);
C = element_prod(A, B); C = element_div(A, B);

transformations


w = conj(u); w = real(u); w = imag(u);
C = trans(A); C = conj(A); C = herm(A); C = real(A); C = imag(A);

Advanced functions

norms


t = norm_inf(v); i = index_norm_inf(v);
t = norm_1(v);   t = norm_2(v); 
t = norm_inf(A); i = index_norm_inf(A);
t = norm_1(A);   t = norm_frobenius(A);

products


axpy_prod(A, u, w, true);  // w = A * u
axpy_prod(A, u, w, false); // w += A * u
axpy_prod(u, A, w, true);  // w = trans(A) * u
axpy_prod(u, A, w, false); // w += trans(A) * u
axpy_prod(A, B, C, true);  // C = A * B
axpy_prod(A, B, C, false); // C += A * B

Note: The last argument (bool init) of axpy_prod is optional. Currently it defaults to true, but this may change in the future. Set the init to true is equivalent to call w.clear() before axpy_prod. Up to now there are some specialisation for compressed matrices that give a large speed up compared to prod.


w = block_prod<matrix_type, 64> (A, u); // w = A * u
w = block_prod<matrix_type, 64> (u, A); // w = trans(A) * u
C = block_prod<matrix_type, 64> (A, B); // w = A * B

Note: The blocksize can be any integer. However, the total speed depends very strong on the combination of blocksize, CPU and compiler. The function block_prod is designed for large dense matrices.

rank-k updates


opb_prod(A, B, C, true);  // C = A * B
opb_prod(A, B, C, false); // C += A * B

Note: The last argument (bool init) of opb_prod is optional. Currently it defaults to true, but this may change in the future. This function may give a speedup if A has less columns than rows, because the product is computed as a sum of outer products.

Submatrices, Subvectors


w = project(u, r); w = project(u, s);
C = project(A, r1, r2); C = project(A, s1, s2);
w = row(A, i); w = column(A, j);

There are to more ways to access some matrix elements as a vector:


matrix_vector_range<matrix_type> (A, r1, r2);
matrix_vector_slice<matrix_type> (A, s1, s2);

Note: These matrix proxies take a sequence of elements of a matrix and allow you to access these as a vector. In particular matrix_vector_slice can do this in a very general way. matrix_vector_range is less useful as the elements must lie along a diagonal.

Example: To access the first two elements of a sub column of a matrix we access the row with a slice with stride 1 and the column with a slice with stride 0 thus:
matrix_vector_slice<matrix_type> (A, slice(0,1,2), slice(0,0,2));

Speed improvements

If you know for sure that the left hand expression and the right hand expression have no common storage, then you can tell ublas that there is no aliasing:


noalias(C) = prod(A, B);

Most often the right hand side of an assignement is constant. So you can give your compiler a hint to use const member function even if A is mutable. (MATRIX is the type of A.) This cast drastically reduces the access time of sparse matrix elements, since no temporary sparse element proxies are created.


C = static_cast<const MATRIX&> A;

Copyright (©) 2000-2004 Joerg Walter, Mathias Koch, Gunter Winkler, Michael Stevens
Permission to copy, use, modify, sell and distribute this document is granted provided this copyright notice appears in all copies. This document is provided ``as is'' without express or implied warranty, and with no claim as to its suitability for any purpose.

Last revised: 2004-07-05