We present a new algorithm to calculate the thermodynamic quantities of three-dimensional (3D) classical statistical systems, based on the ideas of the tensor product state and the density matrix renormalization group. We represent the maximum-eigenvalue eigenstate of the transfer matrix as the product of local tensors that are iteratively optimized by the use of the "vertical density matrix" formed by cutting the system along the transfer direction. This algorithm, which we call vertical density matrix algorithm (VDMA), is successfully applied to the 3D Ising model. Using the Suzuki-Trotter transformation, we can also apply the VDMA to 2D quantum systems, which we demonstrate for the 2D transverse field Ising model.