In recent years, community discovery in heterogeneous networks has gradually become a research hotspot. However, most of the existing methods for discovering non-overlapping or overlapping communities only take one single type of information network into account, and cannot be applied to heterogeneous networks containing multi-mode entities and their multi-dimensional relationships. Presently as a new emerging heterogeneous network, location based social network (LBSN) is attracting more and more attention from social network field. How to effectively discover the hidden complex community structures with multi-dimensional relationships in LBSN, is a very challenging problem for current researchers. Therefore, a community discovery method called multi-relational nonnegative matrix factorization (MRNMF) is proposed that integrates both user and location entities and fuse their multidimensional relationships in LBSN. This method establishes a joint clustering objective function based on nonnegative matrix factorization (NMF), and considers the effect of multi-dimensional factors such as user social relations, user-location check-ins and features of points of interests (POIs). The merits are that not only obtaining accurate user fuzzy communities, but also getting closely related clusters of POIs, which can effectively alleviate data sparse problem in recommendations. The experimental results on two real LBSN datasets show that the proposed method MRNMF has better recommendation performance than other traditional methods in the dual recommendations for POIs and users.