The fractality of complex networks has attracted much attention with extensive investigations over the last 15 years. As a generalization of fractal analysis, multifractal analysis (MFA) is a useful tool to systematically describe the spatial heterogeneity of both theoretical and experimental fractal patterns. One of the widely used methods for fractal analysis is box-covering. It uses the minimum number of covering boxes to calculate the fractal dimension of complex networks, and is known to be NP-hard. More severely, in comparison with fractal analysis algorithms, MFA algorithms have much higher computational complexity. Among various MFA algorithms for complex networks, the sandbox MFA algorithm behaves with the best computational efficiency. However, the existing sandbox algorithm is still computationally expensive. Thus, so far it has only been applied to small-scale complex networks of the size of about tens of thousands of nodes. It becomes challenging to implement the MFA for large-scale networks with tens of millions of nodes. It is also not clear whether or not MFA results can be improved by a largely increased size of a theoretical network. To tackle these challenges, a computationally-efficient sandbox algorithm (CESA) is presented in this paper for MFA of large-scale networks. Distinct from the existing sandbox algorithm that uses the shortest-path distance matrix to obtain the required information for MFA of complex networks, our CESA employs the breadth-first search (BFS) technique to directly search the neighbor nodes of each layer of center nodes, and then to retrieve the required information. Our CESA's input is a sparse data structure derived from the compressed sparse row (CSR) format designed for compressed storage of the adjacency matrix of large-scale network. A theoretical analysis reveals that the CESA reduces the time complexity of the existing sandbox algorithm from cubic to quadratic, and also improves the space complexity from quadratic to linear. MFA experiments are performed for typical complex networks to verify our CESA. The CESA is demonstrated to be effective, efficient and feasible through the MFA results of (u,v)-flower model networks from the 5th to the 12th generations. It enables us to study the multifractality of networks of the size of about 11 million nodes with a normal desktop computer. Furthermore, we have also found that increasing the size of (u,v)-flower model network does improve the accuracy of MFA results. Finally, our CESA is applied to a few typical real-world networks of large scale.