The partition function of the six-vertex model on a square lattice with domain wall boundary conditions (DWBC) is rewritten as a hermitean one-matrix model or a discretized version of it (similar to sums over Young diagrams), depending on the phase. The expression is exact for finite lattice size, which is equal to the size of the corresponding matrix. In the thermodynamic limit, the matrix integral is computed using traditional matrix model techniques, thus providing a complete treatment of the bulk free energy of the six-vertex model with DWBC in the different phases. In particular, in the anti-ferroelectric phase, the bulk free energy and a subdominant correction are given exactly in terms of elliptic theta functions.