
==== Front
J Res Natl Inst Stand Technol
J Res Natl Inst Stand Technol
JRES
Journal of Research of the National Institute of Standards and Technology
1044-677X
2165-7254
[Gaithersburg, MD] : U.S. Dept. of Commerce, National Institute of Standards and Technology

10.6028/jres.125.032
125032
Article
Parallel Generalized Real Symmetric-Definite Eigenvalue Problem
Sims James S. jim.sims@nist.gov
1
Ruiz María Bélen maria.belen.ruiz@fau.de
2
1 National Institute of Standards and Technology, Gaithersburg, MD 20899, USA
2 Friedrich-Alexander-University Erlangen-Nürnberg Egerlandstraße 3, D-91058, Erlangen, Germany
2020
06 11 2020
125 12503222 10 2020
2020
https://creativecommons.org/publicdomain/zero/1.0/ The Journal of Research of the National Institute of Standards and Technology is a publication of the U.S. Government. The papers are in the public domain and are not subject to copyright in the United States. Articles from J Res may contain photographs or illustrations copyrighted by other commercial organizations or individuals that may not be used without obtaining prior approval from the holder of the copyright.
A computationally fast Fortran 90+ quadruple precision portable parallel GRSDEP (generalized real symmetric-definite eigenvalue problem) package suitable for large (80,000 x 80,000 or greater) dense matrices is discussed in this paper.

Key words:

dense matrix eigensolver
high precision
parallel
National Institute of Standards and Technology9999-NIST
==== Body
pmcSoftware DOI: https://doi.org/10.18434/mds2-2293

1 Introduction

Few-electron atomic and molecular systems continue to play an important role as a testbed for understanding fundamental physical and chemical theories. Non-relativistic ground state energies and wave functions for few-electron atoms are primary inputs for computations of atomic properties such as electron affinities and ionization potentials, as well as for computing perturbative corrections due to relativity and quantum electrodynamics [1]. Furthermore, fundamental quantum theory dictates that, in general, uncertainties in wave functions scale as the square root of uncertainties in corresponding energies. Therefore, when using convergence of energies as a metric for convergence of associated wave functions, very high accuracy in the former must be achieved to guarantee accuracy in the later.

To achieve this high accuracy (real(24) arithmetic is necessary for 2 electron systems and real(16) for more than 2 electron systems) requires explicitly correlated wave functions which in turn lead to dense matrices with sizes (N) that are large compared to most dense matrix problems. This in turn leads to linear dependency issues that are solved by going to higher precision, all of which means that calculations need to be parallelized for both computational speed and to spread the memory out across a cluster. Hence parallel packages like ScaLAPACK, which is available only in real(8), cannot be used as they stand and the packages are too big to contemplate conversion to real(16). This has led to the development of the portable parallel GRSDEP (generalized real symmetric-definite eigenvalue problem) package which is discussed in this paper.

2 Parallel GRSDEP

Computational Quantum Chemistry involves solving SchrÖdinger's equation for Ne electrons

Hψ(X1, X2, ..., XNe) = Eψ(X1, X2, ..., XNe), (1)

where Xi = (ri, ξi) is the combined space-spin coordinate for electron i, ψ(X1, X2, ..., XNe) represents a wave function describing a given quantum mechanical system, E is the energy of the system, and H is the Hamiltonian operator specifying components of the energy included in the model.

The wave function given by Equation (1), ψ(X1, X2, ..., XNe), is a linear combination of terms ΦK, where the coefficients cK are those which minimize the total energy, E, given by

E=<ψℋ|ψ><ψ|ψ>=ΣK,LcKcLHKLΣK,LcKcLSKL, (2)

where

HKL =< ΦK |ℋ|ΦL >; SKL =<ΦK|ΦL > . (3)

The condition for the energy to be an extremum, δE = 0, is the well-known matrix eigenvalue (secular) equation:

∑L=1NcLHKL=E∑L=1NcLSKL(4)

Solving this equation is equivalent to solving the N-dimensional generalized eigenvalue problem (GEVP)

Hc = ESc, (5)

where H and S are both symmetric, have matrix elements HKL and SKL given by Equation (3), S is positive definite and N is the number of configurations ΦK in ψ. Equation (5) has N solutions defined by the pairs {(E, c)i, i = 1, N}, each of them fulfilling the variational principle ([2]) for the ith quantum state of the system.

As the Hamiltonian operator H is hermitian, the matrices involved are symmetric (which makes this a generalized real symmetric-definite eigenvalue problem (GRSDEP)), so it is sufficient to compute only the upper or lower triangle of the matrix (an important CPU time and memory saving feature). For use in quantum chemical calculations1 several iterative procedures have been devised for Configuration Interaction(CI) [3-5] and for molecular elliptic coordinate calculations [6]. As pointed out by Lüchow and Kleindienst [7], all of these show slow convergence in the more accurate explicitly correlated calculations and cannot be applied successfully if high precision, spectroscopic accuracy is desired. This slow convergence is probably due to the large off-diagonal matrix elements in the (dense) matrices. For high precision, real(16) arithmetic is necessary for molecules and N > 2 electron atoms, and real(24) arithmetic is necessary for 2 electron atoms. Hence parallel packages like ScaLAPACK, which is available only in real(8), cannot be used as they stand and the packages are too big to contemplate conversion to real(16). For these reasons a portable parallel GRSDEP package was developed.

2.1 Inverse Iteration

Inverse iteration (also called shifted inverse power method in this case) is one of the fastest methods for getting selected roots of the generalized eigenvalue problem, particularly if one has a pretty good idea of what the eigenvalue should be, as is the case in our calculations [8, 9]. When a good starting approximation for an eigenvalue is not available, one uses other methods such as Givens [10] to obtain the starting approximation.

Letting A = H -E0S, inverse iteration, developed first by Wielandt [11], is the power method [12] applied to A-1 = (H -E0S)-1 in the generalized eigenvalue problem, where A-1 is the inverse of A. E0 is some starting approximation for the eigenvalue of interest. If x is an eigenvector of H with associated eigenvalue E (so that Hx = ESx), then subtracting E0S from both sides,

(H -E0S)x = (E -E0)Sx, (6)

x = (E -E0)(H -E0S)-1Sx, (7)

(H -E0S)-1Sx = (E -E0)-1x. (8)

If v0 is not an eigenvector of H it can be expanded in the eigenvectors of H as follows

v0=∑i=1Ncixi(9)

i.e., v1 is the vector obtained by operating with A-1S on v0, v2 is obtained by operating with A-1S on v1, etc.

Letting A = H -E0S, the power method applied to A-1discussed in the Appendix) generates the sequence {vn} recursively, using

vn+1 = A-1Svn. (10)

Equation (10) is the basic iteration formula. Since it is an inverse, the eigenvalue ultimately converged upon will be a LOWER bound to the exact eigenvalue of this state, as required.

The inverse of A is not actually computed (which would not be efficient). Rather, as explained in the next (sub)section, L (and D) are computed and stored in a efficient form for subsequent use in an LT D-1L operation on Svn.

The equation which is used (see Appendix) to monitor the convergence to the final eigenvalue E is

E=E0+vn+1·Svnvn+1·Svn+1(1)

The convergence of vn to x is faster the closer the trial E0 is to E. Hence it is desirable to have a pretty good idea of what the eigenvalue is. It is generally sufficient to know the eigenvalue to 5 digits. When that is not the case, one can use NAG [13] or LAPACK [14] libraries with smaller expansions to get the starting value.2

2.2 Strategy for using Inverse Iteration

The (efficient) strategy for the inverse iteration without actually calculating A-1 = (H -E0S)-1 is as follows.

By construction a lower triangular matrix L such that LA = U, an upper triangular matrix is found. Right multiply by LT , where LT is the transpose of L, to get

LALT = ULT = LUT = D, (12)

where D is a diagonal matrix since LALT is symmetric. Here the facts that the set of upper (lower) triangular matrices is closed under multiplication and that the intersection of the set of upper triangular matrices and lower triangular matrices is the set of diagonal matrices are used. Now solve for A to get

A = L-1DL-T (13)

and

A-1 = LT D-1L (14)

where we have used the fact that the inverse of the transpose is the transpose of the inverse. Instead of A-1 in Equation (10) LT D-1L is used.

It remains only to show how to determine L and D in the transformation of A to diagonal form. L is assembled by construction, namely,

L = LNLN-1...L3L2 (15)

and

A2 = L2A (16)

A3 = L3A2 (17)

... (18)

Ak = LkAk-1. (19)

L2 is the lower triangular matrix that sets the off-diagonal lower triangular elements of the second row of A to zero, L3 operating on A2 does a similar thing for the third row of A2, and so on. In general, Lk is a matrix with 1s along the diagonal and zeros everywhere else except the lower triangular off-diagonal elements of row k, which are given by

- Ak-1(i, k)/Ak-1(i, i), i = 1, 2, ..., k - 1; (20)

that is, the elements of the k-th column of Ak-1 divided by the corresponding diagonal elements.

Once Ak is obtained column k (above the diagonal) is no longer needed and the off-diagonal elements of row k of Lk can be safely stored there. Once AN is obtained, the diagonal elements constitute the diagonal matrix D. For convenience D is then replaced by D-1. Vector A is now what we call the implicit representation of L. There is no need to multiply out the Ls (an Order N3 operation); each operation on x by an Lk is at most an order N operation.

Note that the basic operation (Ak = LkAk-1) in forming the Aks is the Fortran 90+ DOT PRODUCT intrinsic function, which hopefully will be implemented efficiently. L = LNLN-1...L3L2 also involves only the DOT PRODUCT. LT y (y = D-1L) on the other hand uses the QAXPY BLAS operation, which also can be efficiently implemented for portability and extended precision.

A is not stored as an N by N array but rather as a vector of upper triangular columns, which are efficiently accessed by the DOT PRODUCT function (and similarly for QAXPY).

One thing to note is that with the large matrices that one deals with in parallel problems, memory access patterns can be extremely important. The initial coding got very slow as large N was approached. An essentially trivial reordering of the operations involved with the Ls fixed the problem (a ROW form of L and a COLUMN form for H and S are now used) and led to substantial speedups for the large N case. All of which emphasizes that the construction of L is by no means unique and that this one is a very good one.

2.3 Parallelization

To parallelize one simply maps A to the several processors. A cyclic range is used to distribute elements of each row of the matrix A over the processors in a cluster. Each column is stored in its entirety on a particular processor. A cyclic "round robin" assignment of columns to the processors is made, which yields an approximately equal distribution of matrix elements to processors. This insures good load balancing (assuming homogeneous processors) and produces better load balancing than a block assignment (as in ScaLAPACK [16]) would.

In the parallel version of GRSDEP, MPI [17] is used to run the same program on multiple processors (on the same or different hosts) and each processor generates H and S only for the columns of H and S it owns.

The steps of the (parallel) inverse iteration are

1. Calculate A = H -ES.

On each processor, calculate the As for the columns of the matrices that belong to this processor.

2. Calculate Sv0 using the vector 1=(1,1,...,1) as the starting vector v0.

On each processor, calculate the partial vectors bi = Si · 1 (a row sum, 1 is a vector of all ones) and send them to the root process (process 0). Process 0 then sums them up to get b = S · 1(= Sv0).

3. Find the transformation of A to the implicit representation of L.

Calculate the partial Ls of A with the partial Ds stored in the diagonal elements. Set k = 0.

4. Begin iteration.

(a) If (n ≠ 0) Calculate bn = Svn.

(b) Solve vn+1=LT D-1Lbn (Call Ltdm1L to do LT D-1Lbn, returning the result as x which is vn+1).

5. Calculate E = E0 + vn+1·SvnSvn+1·vn+1.

6. Stop condition fulfilled?

(E changed negligibly between last two iterations?)

NO: set vn+1 = vn+1rwhere r =vn+1·Svn+1; n ← n + 1; go to 4 (a).

YES: stop.

One other point to discuss is how to handle the Ls and Ds. Do we keep them distributed during the operation with L on x, or are the distributed Ls gathered into processor 0 and the L operations carried out on a single processor (that is serially; but even with multiple processors the operation by L is still strictly serial). The problem with keeping L on a single processor is that it unbalances the memory requirements. So L (on top of A) is kept distributed and the vectors are passed to the processors in a round-robin fashion. So for the Lx operation

1. process 0 applies its Ls to x and passes the resulting vector to process 1.

2. process 1 applies its Ls to the x from process 0, then passes result to process 2.

3. process 2 applies its Ls to the vector received from process 0 and so on.

Except for root (P0) a processor Pi receives an x from Pi-1 and sends the x it calculates to Pi+1, except if i + 1 is greater than last, in which case x will be sent to P0, which will already have issued a receive. All processors except P0 start out with a receive from processor myrank - 1 and wait until a message arrives. For LTx one uses a similar logic except one starts with Plast and works down to P0.

3 Summary of the Timing Tests

This section is an attempt to measure the effectiveness of the parallelization of the GRSDEP by calculating the time it takes (in seconds (s)) to compute the energy of a 4190 term H2 wave function with a resultant total energy of -1.1744 7571 hartrees (Ha) [9] on differing numbers of processors. Table 1 gives the times (in seconds (s)). All results were obtained using real(16) precision (quadruple precision or QP, 128-bit, ≈ 32 digits) floating point arithmetic on the National Institute of Standards and Technology (NIST)'s 16908 processor Linux cluster. The Message Passing Interface (MPI) Standard [17] was used to parallelize the code.

Table 1 The time it takes to compute the ground state energy of the hydrogen molecule as a function of the number of processors utilized on a Linux cluster.

Number of
Processors	LD
time (s)	Ltdm1L
time (s)	Total
time (s)	
1	1392	26	1428	
2	697	26	728	
4	351	26	380	
8	177	26	205	
16	89	26	117	
32	46	26	73	

The LD step in the table refers to the time it takes to find L and D and gives a speedup of 30 on 32 processors, which is excellent! The 1 to 32 run speeds up by a factor of 20, a 23 minute job completing in a little over a minute. While performance can depend on the amount of memory available, the size of the problem (N), how one calculates LT, how elements of A are stored, etc., the speedup is ultimately limited by the two sequential steps, calculating the matrix elements and the Ltdm1L (LT D-1L) step.

4 Conclusions

While the parallel GRSDEP has proven quite useful, having been utilized for investigating the following systems:

⏺ 2-electron molecule: ground state of H2,

⏺ 2-,3-, 4-, and 5-electron atoms and their ions,

and producing some of the lowest (i.e. most accurate) variational energies ever reported, the lack of pivoting in this solver needs to be discussed.

The solver (shifted inverse iteration) used is one of the fastest methods for getting selected roots of the generalized eigenvalue problem. The only computational problem is the computation of the implicit representation of L. The addition of partial and/or complete pivoting in the use of inverse iteration complicates the code as the matrices go from triangular to square and much of the speed advantage of inverse iteration is lost. Equation (10) is also

Avn+1 = Svn (21)

which is the familiar linear system of the form Ax = b for which LAPACK has specific routines with pivoting. One of these which became available in version 3.5 utilizes the bounded Bunch-Kaufman (rook) pivoting algorithm [18] which is more accurate than partial pivoting, with comparable costs. Since LAPACK is available in Fortran, it was not hard to produce real(16) versions of those routines in place of step 3 and 4.b in the inverse iteration algorithm (sequential case) outlined in Section 2.3 above, and hence test the solver-produced results reported in, for example [19]. It should be noted that the Ritz estimate

E=E0+vn+1·Hvn+1vn+1·Svn+1. (22)

is not used to monitor convergence to the final eigenvector E, rather it is Equation (11), derived in the Appendix, which is used in step 2.3.5. This has the advantage that the original H matrix does not need to be saved, freeing that storage space up to be overwritten by L, which then can be used to solve for a principal minor of dimension N' with very little additional CPU time.

In all cases the results reported in [19] were reproduced to all digits reported in that work. However it took 9.5 days to reproduce the result for the largest (N = 39,381) matrix (which took only 176 minutes on 120 processors), obviously not fast enough for a production code even for matrices of this size or smaller, much less for the sizes computed (up to N ≈ 80 K) in earlier work [20]. Public domain codes are available for inverse iteration with partial and complete pivoting, but they do not meet the need for a parallel code which is fast, efficient, numerically stable, and of high enough precision (real(16)) to enable calculations which go significantly beyond the currently investigated sizes. A real(16) parallel version of the Rook pivoting algorithm, such as the one(s) of Strazdins [21], could greatly facilitate this endeavor.

5 Software Specifications

NIST Operating Unit	NIST Information Technology, Applied and Computational Mathematics,
Scientific Applications and Visualization	
Category	parallel GEVP (Generalized Eigenvalue Problem)	
Targeted Users	High precision atomic and molecular science	
Operating Systems	All	
Programming Language	Fortran 90+	
Inputs/Outputs	See the HowTo included with the software: https://doi.org/10.18434/mds2-2293	
Documentation	Included with the software: https://doi.org/10.18434/mds2-2293	
Disclaimer	https://www.nist.gov/director/licensing	

Acknowledgments

This work would not have been possible without Stanley Hagstrom, who was always willing to share his seemingly infinite knowledge of chemistry and computer science with us. One of us would like to thank William George for MPI Support, Carl Spangler and Denis Lehane for parallel systems support, and Judy Devaney Terrill and Jules Becker for encouraging the publication of this work.

About the authors: James S. Sims, a computational scientist, is a former staff member (retired) and now a Research Associate in the Scientific Applications and Visualization Group, Applied and Computational Mathematics Division of the NIST Information Technology Laboratory.

María Belén Ruiz is a Privatdozentin (independent lecturer and researcher) in the Department of Chemistry and Pharmacy at the Friedrich-Alexander-University Erlangen-NuA"rnberg in Erlangen, Germany.

The National Institute of Standards and Technology is an agency of the U.S. Department of Commerce.

6 Appendix: The power method applied to A = H -ES

Letting v1 be the vector which results from applying A-1S to v0,

v1 = (H -E0S)-1Sv0 = (H -E0S)-1∑i=1Ncixi(23)

= ∑i=1Nci(H -E0S)-1Sxi = ∑i=1Nci (Ei -E0)-1xi. (24)

Simliarly, let v2 be the vector which results from applying A-1Sv to v1, and similarly for v3, etc.

v2 = (H -E0S)-1Sv1 = ∑i=1Nci (Ei -E0)-2xi (25)

...

vn = (H- E0S)-1Svn -1 = ∑i=1Nci (Ei -E0)-nxi. (26)

From Equation (26) it is clear that if (Ek -E0) is small, after a few iterations the right hand side will be proportional to the eigenvector xk (all terms but the i = k term will be small), so that

vn ≈ ck(Ek -E0)-nxk. (27)

Setting Ek = E and applying A-1S again to vn, one gets

vn+1 = A-1Svn, (28)

vn+1 = ck(E -E0)-(n+1)xk, (29)

vn+1 = (E -E0)-1ck(E -E0)-nxk, (30)

vn+1 = (E -E0)-1vn, (31)

(E -E0)vn+1 = vn. (32)

Applying S to both sides and then taking the dot product of both sides with vn+1,

(E -E0)vn+1 · Svn+1 = vn+1 · Svn,(33)

E-E0=vn+1·Svnvn+1·Svn+1 (34)

E=E0+vn+1·Svnvn+1·Svn+1. (35)

Which is Equation (11).

1 With applicability in other fields, such as atomic, molecular and optical physics, as well.

2 When we have no "guess" for the desired eigenvalue, we use a truncated wave function and solve the eigenvalue problem by reduction of S to a unit matrix followed by solving the transformed H for the lowest few eigenvalues using Givens [15].
==== Refs
7. References

1 King FW (1999) High-precision calculations for the ground and excited states of the lithium atom. Advances In Atomic, Molecular, and Optical Physics 40 :57-112.
2 Shull H, Löwdin PO (1958) Variational Theorem for Excited States. Physical Review 110 (6 ):1467.
3 Nesbet RK (1965) Algorithms for Diagonalization of Large Matrices. Journal of Chemical Physics 43 :311-312.
4 Davidson ER (1975) The iterative calculation of a few of the lowest eigenvalues and corresponding eigenvectors of large real-symmetric matrices. Journal of Computational Physics 17 :87-94.
5 Shavitt I, Bender CF, Pipano A, Hosteny RP (1973) The iterative calculation of a few of the lowest eigenvalues and corresponding eigenvectors of large real-symmetric matrices. Journal of Computational Physics 11 :90-108.
6 Cheung LM, Bishop DM (1977) The group-coordinate relaxation method for solving the generalized eigenvalue problem for large real-symmetric matrices. Computer Physics Communications 13 :247-250.
7 Lüchow A, Kleindienst H (1993) Stable and efficient algorithm for selected eigenvalues and eigenvectors of the general symmetric eigenproblem. Computers & Chemistry 17 (1 ):61-66.
8 Sims JS, Padhy B, Ruiz MB (2020) Exponentially correlated Hylleraas-configuration interaction (E-Hy-CI) nonrelativistic energy for the 1S ground state of the helium atom. International Journal of Quantum Chemistry 12 October 2020:e26470.
9 Sims JS, Hagstrom SA (2006) High precision variational calculations for the Born-Oppenheimer energies of the ground state of the hydrogen molecule. Journal of Chemical Physics 124 :094101.
10 Ralston and Wilf (1967) Mathematics for digital computers (John Wiley and Sons) Vol. 2, p 94ff.
11 Wielandt H (1944) Ber B4/j/37 Aerodyn Vers-Anst .
12 Ralston A, Rabinowitz P (1978) A First Course in Numerical Analysis (McGraw-Hill, New York).
13 The Numerical Algorithms Group (2020) The NAG Library, Oxford, United Kingdom.
14 Anderson E, Bai Z, Bischof C, Blackford S, Demmel J, Dongarra J, Du Croz J, Greenbaum A, Hammarling S, McKenney A, Sorensen D (1999) LAPACK Users' Guide (Society for Industrial and Applied Mathematics, Philadelphia, PA) 3rd Ed.
15 Moler C, Spangler D (1979) Eispack-based substitute for QCPE subroutine Givens. The National Resource for Computation in Chemistry.
16 Blackford LS, Choi J, Cleary A, D'Azevedo E, Demmel J, Dhillon I, Dongarra J, Hammarling S, Henry G, Petitet A, Stanley K, Walker D, Whaley RC (1997) ScaLAPACK Users' Guide (Society for Industrial and Applied Mathematics, Philadelphia, PA).
17 Message Passing Interface Forum (1994) MPI: A Message-Passing Interface Standard. The International Journal of Supercomputer Applications 8 (3/4 ):159-416.
18 Khabou A, Demmel JW, Grigori L, Gu M (2013) LU Factorization with Panel Rank Revealing Pivoting and Its Communication Avoiding Version. SIAM Journal on Matrix Analysis and Applications 34 (3 ):1401-1429.
19 Sims JS (2020) Hylleraas-configuration interaction (Hy-CI) non-relativistic energies for the 3 1S, 4 1S, 5 1S, 6 1S, and 7 1S excited states of the beryllium atom. Journal of Research of the National Institute of Standards and Technology 125 :125006.
20 Sims JS, Hagstrom SA (2014) Hylleraas-configuration-interaction nonrelativistic energies for the 1S ground states of the beryllium isoelectronic sequence. Journal of Chemical Physics 140 :224312. 24929393
21 Strazdins P (2000) Accelerated methods for performing the LDLT decomposition. ANZIAM 42 :C1328-C1355.
