John E. Pask

Eugene Vecharynski, Chao Yang, John E. Pask

We present an iterative algorithm for computing an invariant subspace associated with the algebraically smallest eigenvalues of a large sparse or structured Hermitian matrix A. We are interested in the case in which the dimension of the invariant subspace is large (e.g., over several hundreds or thousands) even though it may still be small relative to the dimension of A. These problems arise from, for example, density functional theory based electronic structure calculations for complex materials. The key feature of our algorithm is that it performs fewer Rayleigh--Ritz calculations compared to existing algorithms such as the locally optimal precondition conjugate gradient or the Davidson algorithm. It is a block algorithm, hence can take advantage of efficient BLAS3 operations and be implemented with multiple levels of concurrency. We discuss a number of practical issues that must be addressed in order to implement the algorithm efficiently on a high performance computer.

Phanish Suryanarayana, Phanisri P. Pratapa, John E. Pask

We present the Alternating Anderson-Richardson (AAR) method: an efficient and scalable alternative to preconditioned Krylov solvers for the solution of large, sparse linear systems on high performance computing platforms. Specifically, we generalize the recently proposed Alternating Anderson-Jacobi (AAJ) method (Pratapa et al., J. Comput. Phys. (2016), 306, 43--54) to include preconditioning, discuss efficient parallel implementation, and provide serial MATLAB and parallel C/C++ implementations. In serial applications to nonsymmetric systems, we find that AAR is comparably robust to GMRES, using the same preconditioning, while often outperforming it in time to solution; and find AAR to be more robust than Bi-CGSTAB for the problems considered. In parallel applications to the Helmholtz and Poisson equations, we find that AAR shows superior strong and weak scaling to GMRES, Bi-CGSTAB, and Conjugate Gradient (CG) methods, using the same preconditioning, with consistently shorter times to solution at larger processor counts. Finally, in massively parallel applications to the Poisson equation, on up to 110,592 processors, we find that AAR shows superior strong and weak scaling to CG, with shorter minimum time to solution. We thus find that AAR offers a robust and efficient alternative to current state-of-the-art solvers, with increasing advantages as the number of processors grows.

Yunfeng Cai, Zhaojun Bai, John E. Pask et al.

By extending the classical analysis techniques due to Samokish, Faddeev and Faddeeva, and Longsine and McCormick among others, we prove the convergence of preconditioned steepest descent with implicit deflation (PSD-id) method for solving Hermitian-definite generalized eigenvalue problems. Furthermore, we derive a nonasymptotic estimate of the rate of convergence of the \psdid method. We show that with the proper choice of the shift, the indefinite shift-and-invert preconditioner is a locally accelerated preconditioner, and is asymptotically optimal that leads to superlinear convergence. Numerical examples are presented to verify the theoretical results on the convergence behavior of the \psdid method for solving ill-conditioned Hermitian-definite generalized eigenvalue problems arising from electronic structure calculations. While rigorous and full-scale convergence proofs of preconditioned block steepest descent methods in practical use still largely eludes us, we believe the theoretical results presented in this paper sheds light on an improved understanding of the convergence behavior of these block methods.

Gaigong Zhang, Lin Lin, Wei Hu et al.

Recently, we have proposed the adaptive local basis set for electronic structure calculations based on Kohn-Sham density functional theory in a pseudopotential framework. The adaptive local basis set is efficient and systematically improvable for total energy calculations. In this paper, we present the calculation of atomic forces, which can be used for a range of applications such as geometry optimization and molecular dynamics simulation. We demonstrate that, under mild assumptions, the computation of atomic forces can scale nearly linearly with the number of atoms in the system using the adaptive local basis set. We quantify the accuracy of the Hellmann-Feynman forces for a range of physical systems, benchmarked against converged planewave calculations, and find that the adaptive local basis set is efficient for both force and energy calculations, requiring at most a few tens of basis functions per atom to attain accuracy required in practice. Since the adaptive local basis set has implicit dependence on atomic positions, Pulay forces are in general nonzero. However, we find that the Pulay force is numerically small and systematically decreasing with increasing basis completeness, so that the Hellmann-Feynman force is sufficient for basis sizes of a few tens of basis functions per atom. We verify the accuracy of the computed forces in static calculations of quasi-1D and 3D disordered Si systems, vibration calculation of a quasi-1D Si system, and molecular dynamics calculations of H$_2$ and liquid Al-Si alloy systems, where we find excellent agreement with independent benchmark results in literature.