  Within the parameter space of equation of state (EOS) of dense neutron-rich
matter limited by existing constraints mainly from terrestrial nuclear
experiments, we investigate how the neutron star maximum mass
$M_{\rm{max}}>2.01\pm0.04$ M$_\odot$, radius $10.62<R_{\rm{1.4}}< 12.83$ km and
tidal deformability $\Lambda_{1.4}\leq800$ of canonical neutron stars all
together constrain the EOS of dense neutron-rich nucleonic matter. While the
3-D parameter space of $K_{\rm{sym}}$ (curvature of nuclear symmetry energy),
$J_{\rm{sym}}$ and $J_0$ (skewness of the symmetry energy and EOS of symmetric
nuclear matter, respectively) are narrowed down significantly by the
observational constraints, more data are needed to pin down the individual
values of $K_{\rm{sym}}$, $J_{\rm{sym}}$ and $J_0$ with quantified
uncertainties. The $J_0$ largely controls the maximum mass of neutron stars.
While the EOS with $J_0=0$ is sufficiently stiff to support neutron stars as
massive as 2.37 M$_{\odot}$, to support the hyperthetical ones as massive as
2.74 M$_{\odot}$ (composite mass of GW170817) requires $J_0$ to be larger than
its currently known maximum value of about 400 MeV and beyond the causality
limit. The upper limit on the tidal deformability of $\Lambda_{1.4}=800$ from
the recent observation of GW170817 is found to provide upper limits on some EOS
parameters consistent with but far less restrictive than the existing
constraints of other observables studied.
