Scientific Computing Lecture 5 - Matrix Operations 1

Definitions
Addition
Norms
Matrix vector product
Naïve matrix product
Alternative implementations
Blocking

Vector can be considered as m×1 matrix
Matlab does this (everything is a matrix in matlab)
Scalar values are 1×1 matrices
However, vectors are special cases of
matrices concerning representation
Mathematically usually different meaning
Submatrices gained by deleting rows and
columns of matrix

Splitting matrix into submatrices
Partitioned matrix sometimes called block
matrix
In most cases, submatrices of block matrix have
same size

An n×n matrix is called
Diagonal if a_ij = 0 for all i ≠ j
Upper triangular if a_ij = 0 for all i > j
Lower triangular if a_ij = 0 for all i < j
A partitioned matrix is called
Upper block triangular matrix if A_ij = 0 for all i > j

The transposed A^T of an m×n matrix is
n×m matrix where a_ij and a_ij are swapped
Correspondingly the transposed conjugate matrix A^H

A is called symmetric if A=A^T
n A is called hermitian if A=A^H

Matrices are added element-wise
Dimensions must be equal

Is this operation faster or slower then vector
addition?
Is it better to iterate first over the
rows or first over the columns?

Adding matrices of dimensions m×n takes
2 m×n loads
m×n stores
m×n additions
Recall that the additions can be in 1 or 2
processor cycles
Loading and storing is more expansive,
especially from main memory
But also from L2 (or L3 if exist)

Adding vectors of length l = m×n takes
2 l loads
l stores
l additions
Thus same arithmetic and memory operations are
performed
Should run at the same speed, unless …
See next slides

In C and C++ 2D arrays are stored row major
All elements of one row are stored together

In Fortran 2D arrays are stored column major
All elements of one column are stored together
If you define your own matrix type with your own
indexing you use either one
What is done in some libraries, like ???

Loop nesting should fit storing scheme
For row major matrices iterating first over rows
For column major first over columns
Stride-1 memory access
Why is this important?

Matrices are operators on vectors
Matrix vector product Ax can be considered as
applying the operator A on x
Therefore most matrix norms are operator
norms
n Definitions of vector norms be applied to the
image of the operation

||.||_p for p = 1, 2, and ∞ are the most
interesting
p = 2 difficult to compute
p = 1 and ∞ much simpler

The first is the maximum of column sums
The latter is the maximal row sum

Number of columns must be equal to vector
dimension
Result is a vector
Dimension is number of rows

i-th element of result is dot product of i-th
row and multiplied vector

Dot product is loop over matrix rows
What corresponds to a loop over columns of
the matrix?
Scaled vector addition

Is this faster or slower than dot product?

Operand vector: n loads
Result vector: m stores
Matrix: mn loads
Multiplications: mn
Additions: m(n-1)
Which is the most expensive?

Number of additions, multiplications and loads for
the matrix are almost equal
Loading from memory is much slower then basic
arithmetic
Assume we have a row-major matrix and a cache
line of s_c double
Loads are composes of:
Vector addition: mn memory loads
Dot product: mn memory loads and mn loads
from cache
Safe to assume
For column-major matrices it is the opposite