We examine the problem of completely positive (CP) factorization for a given completely positive matrix. We are inspired by the idea of Groetzner and Dür in 2020, wherein the CP factorization problem can be reformulated as an equivalent feasibility problem containing an orthogonality constraint. We begin by solving this feasibility problem through the use of a special case of the Riemannian smoothing steepest descent method proposed by Zhang et al. in 2021. To apply it to the CP factorization problem, we employ a smooth approximation function, named LogSumExp, as the maximum function. Numerical experiments show the efficiency of our method especially for large-scale factorizations.