  The three-dimensional intrinsic shape of a galaxy and the mass of the central
supermassive black hole provide key insight into the galaxy's growth history
over cosmic time. Standard assumptions of a spherical or axisymmetric shape can
be simplistic and can bias the black hole mass inferred from the motions of
stars within a galaxy. Here we present spatially-resolved stellar kinematics of
M87 over a two-dimensional $250\mbox{$^{\prime\prime}$} \times
300\mbox{$^{\prime\prime}$}$ contiguous field covering a radial range of 50
pc-12 kpc from integral-field spectroscopic observations at the Keck II
Telescope. From about 5 kpc and outward, we detect a prominent 25
$\mathrm{km~s}^{-1}$ rotational pattern, in which the kinematic axis
(connecting the maximal receding and approaching velocities) is $40^\circ$
misaligned with the photometric major axis of M87. The rotational amplitude and
misalignment angle both decrease in the inner 5 kpc. Such misaligned and
twisted velocity fields are a hallmark of triaxiality, indicating that M87 is
not an axisymmetrically shaped galaxy. Triaxial Schwarzschild orbit modeling
with more than 4000 observational constraints enabled us to determine
simultaneously the shape and mass parameters. The models incorporate a radially
declining profile for the stellar mass-to-light ratio suggested by stellar
population studies. We find that M87 is strongly triaxial, with ratios of
$p=0.845$ for the middle-to-long principal axes and $q=0.722$ for the
short-to-long principal axes, and determine the black hole mass to be
$(5.37^{+0.37}_{-0.25}\pm 0.22)\times 10^9 M_\odot$, where the second error
indicates the systematic uncertainty associated with the distance to M87.
