P573 Scientific Computing Lecture 8 Sparse Matrices

1
P573Scientific ComputingLecture 8 - Sparse
Matrices

• Peter Gottschling
• pgottsch_at_cs.indiana.edu
• www.osl.iu.edu/pgottsch/courses/p573-06

2
Overview

• Structured sparse matrices
• Band matrices
• Matrix vector product
• Unstructured matrices
• Matrix vector product
• Other operations
• Insertion

3
Tridiagonal Matrices

• Diagonal and off-diagonals are only non-zero
• For instance 1D-Poisson equation
• Can be stored as vectors

4
Band Matrices

• Extension of tridiagonal matrices
• An (l, u)-m?n band matrices has
• Only non-zeros in l diagonals below the main
  diagonal
• Only non-zeros in u diagonals above the main
  diagonal
• Formally

5
Matrix of 2D-Poisson Equation

• Off-diagonals with non-zeros only doesnt need to
  stored
• Instead of 0 diagonals we store offsets of
  non-zero off-diagonals
• In example of 3x3 domain -3, -1, 0, 1, 3
  o0, ..., on

6
Matrix Vector Multiplication yAx

• struct_matrixmult (const vector x, vector y)
• for (int i 0 i lt max(-off0, 0) i)
• for (int j 0 j lt n j)
• if (ioffj gt 0) yi aij
  xioffj
• for (int i max(-off0, 0) i lt min(N,
  N-offn-1) i)
• for (int j 0 j lt n j)
• yi aij xioffj
• for (int i min (N, N-offn-1).)
• Why can ifs be expensive?

7
Matrix-free Matrix Vector Multiplication

• For 9?9-Matrix of Poisson example
• We know the matrix already why we need to store
  it?
• sample_matrixmult (const vector x, vector y)
  
• for (int i 0 i lt 9 i)
• yi 4.0xi
• if (i gt 2) yi -1.0xi-3 //
  northern neig. 2D
• if (i lt 7) yi - xi3 //
  southern neig. in 2D
• if (i 3) yi - xi-1 //
  western neig. in 2D
• if ((i2) 3) yi - xi1 //
  eastern neig in 2D

8
Unstructured Matrices

• If there is no structure or to hard to explore
• For all non-zeros, row, column and value need to
  be stored
• Somehow
• Example 77-Matrix in picture (0,0,1), (0,5,2),
  (1,1,3), ...
• Which options we have to store this information?
• Can we save some memory?

9
Coordinate Matrix

• Use one vector for each elements row, column and
  value
• Vector length is number of non-zeros
• Which advantages and disadvantages has this
  format?

10
CRS Format

• Compressed Row Storage (sometimes CSR)
• rowi first non-zero in row i (0...6N-1)
• rowN nnz past-end index
• columnk column of k-th non-zero
• aijk value of k-th non-zero

11
yAx with Compressed Rows

• crs_matrixmult (const vector x, vector y)
• for (int i 0 i lt N i)
• yi 0
• for (int n rowi, next rowi1 n lt
  next n)
• yi aijn xcolumnn

12
Possible Performance Improvements

• Unroll outer loop
• More independent operations
• Unrolling inner loop only useful if enough
  non-zero per row
• Temporary variable avoids cache invalidation
• crs_matrixmult (const vector x, vector y)
• for (int i 0 i lt N i)
• temp 0
• for (int n rowi, next rowi1 n lt
  next n)
• temp aijn xcolumnn
• yi temp

13
CCS Format

• Compressed Column Storage
• Sometimes CSC
• Roles of rows and columns are interchanged
• In MTL same implementation
• How would you implement a matrix vector product
  with this format?

14
Matrix Product vs. Matrix Vector Product

• Assume we have to multiply matrices A and B
• And some vectors, e.g. y 2ABx
• Which program fragment is more efficient?
• sparse_matrix C mult(A, B, C)
• for (int i 0 i lt 10000 i)
• C.mult(x, z) y 2z
• for (int i 0 i lt 10000 i)
• A.mult(x, z) B.mult(z, w) y 2w

15
Symmetric Matrices

• Store only upper triangular
• Lower is mirrored
• Adapt matrix vector product
• for (int i 0 i lt N i)
• yi 0
• for (int n rowi, next rowi1 n lt
  next n)
• yi aijn xcolumnn
• if (i ! columnn) ycolumnn aijn
  xi
• Can we traverse easily all elements in a row?
• Or in a column?
• For instance to implement matrix addition
• Same applies for anti-symmetric or Hermitian
  matrices
• How to generalize the implementation?

16
Summary

• Different representations
• If matrix has structure, exploit it for
  performance
• For unstructured matrices CRS format favorable
• Sparse matrix multiplication is rare
• Replacing by sequence of matrix vector product in
  most cases better
• Most important operation for sparse matrices is
  matrix vector product
• One of the most important components in
  simulations
• If not the most important performance-wise