We report an implementation of the LDA+U method based on the state-of-the-art linear combination of pseudo-atomic orbital (LCPAO) method, which is suitable for large-scale O(N) electronic structure calculations based on the density functional theory. By introducing a dual representation of the occupation number matrix instead of the on-site or full representations, the LDA+U formalism is refined to be consistent with a nonorthogonal LCPAO basis in regard to the sum rule of the total number of electrons. For typical transition metal oxide bulk systems, the band gap, magnetic moment, and detailed electronic structures are investigated with the different choices of basis orbitals and effective U values as well as the definition of the occupation number matrix. The results are in good agreement with previous theoretical and experimental studies, indicating that the proposed LDA+U scheme combined with the O(N) method is a quite promising approach for the study of large-scale correlated material systems consisting of localized electrons. We discuss the electronic structure and magnetic properties of (NiO)(m)/(CoO)(n) superlattices as an application of our method.