This is a generalized version of bivariate P-splines for covariance smoothing of sparse functional or longitudinal data. It uses tensor product B-spline basis functions and employs a differencing penalty on the associated parameter matrix. The only smoothing parameter in the method is selected by leave-one-subject-out cross-validation and is implemented with a fast algorithm.