(Click on title slide to access entire slide deck.)

Watch the lecture here.

Lecture 6

In this lecture we develop a simple model for CPU operation that we use to motivate a sequence of optimizations for matrix-matrix product.

The matrix class we developed in Lecture 5 is straightforward:

#include <vector>

class Matrix {
public:
  Matrix(size_t M, size_t N) : num_rows_(M), num_cols_(N), storage_(num_rows_ * num_cols_) {}

        double& operator()(size_t i, size_t j)       { return storage_[i * num_cols_ + j]; }
  const double& operator()(size_t i, size_t j) const { return storage_[i * num_cols_ + j]; }

  size_t num_rows() const { return num_rows_; }
  size_t num_cols() const { return num_cols_; }

private:
  size_t              num_rows_, num_cols_;
  std::vector<double> storage_;
};

But, like the vector class we developed in Lectures 4 and 5, this class is extremely efficient and useful.

One of the main design principles behind our matrix and vector classes is that we want to be able to implement algorithms for them that only use their external interfaces – i.e., operator(). Algorithms developed this way do not depend on any internal implementation details of the matrix and vector classes.

For example, a matrix-vector product algorithm might look like the following:

void   matvec(const Matrix& A, const Vector& x, Vector& y) {
  for (size_t i = 0; i < A.num_rows(); ++i) {
    for (size_t j = 0; j < A.num_cols(); ++j) {
      y(i) += A(i, j) * x(j);
    }
  }
}

Note that this computes ( y \leftarrow y + A \times x ). This is actually slightly more useful (being more general) as we can compute ( y \leftarrow A \times x ) simply by passing in ( y = 0 ).

Operator syntax for computing ( y \leftarrow A \times x ) is then

Vector operator*(const Matrix& A, const Vector& x) {
  Vector y(A.num_rows());
  matvec(A, x, y);
  return y;
}