Lecture 6: CPUs¶

(Originally recorded 2019-04-18)

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 well as being a common oepration in numerical linear algebra. We can still 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;
}

Note also that the matvec() function does not return anything – it stores its results into one of its arguments. The operator*() function, on the other hand, does create a new Vector and returns that.