Quasi-Monte Carlo (QMC) rules 1/N Sigma(N-1)(n=0) f(y(n)A) can be used to approximate integrals of the form integral([0,1]s) f(yA) dy, where A is a matrix and y is a row vector. This type of integral arises, for example, from the simulation of a normal distribution with a general covariance matrix, from the approximation of the expectation value of solutions of PDEs with random coefficients, or from applications from statistics. In this paper we design QMC quadrature points y(0), . . ., y(N-1) is an element of [0, 1](s) such that for the matrix Y = (y(0)(inverted perpendicular), . . ., y(N-1)(inverted perpendicular))(inverted perpendicular) whose rows are the quadrature points, one can use the fast Fourier transform to compute the matrix-vector product Y a(inverted perpendicular), a is an element of R-s, in O(N log N) operations and at most 2(s - 1) extra additions. The proposed method can be applied to lattice rules, polynomial lattice rules, and a certain type of Korobov p-set and even works if the point set y(0), . . ., y(N-1) is transformed to another domain U subset of R-s by a coordinatewise mapping phi which is the same in each coordinate. The approach is illustrated computationally by three numerical experiments. The first test considers the generation of points with normal distribution and a general covariance matrix, the second test applies QMC to high-dimensional, affine-parametric, elliptic PDEs with uniformly distributed random coefficients, and the third test addresses finite element discretizations of elliptic PDEs with high-dimensional, log-normal random input data. All numerical tests show a significant speedup of the computation times of the fast QMC matrix method compared to a conventional implementation as the dimension becomes large.