The Yoshimine Sorting Algorithm


T. Daniel Crawford


Center for Computational Quantum Chemsitry, Department of Chemistry,
The University of Georgia, Athens, Georgia 30602-2525
and
Institute for Theoretical Chemistry, Departments of Chemistry and Biochemistry,
The University of Texas, Austin, Texas 78712-1167

Quantum chemical methods generally require storage and manipulation of large quantities of data such as two-electron integrals, configuration interaction coefficients, cluster amplitudes, etc. Often, the algorithms designed to make use of this data demand that it be ordered in a particular manner. For example, in coupled-cluster theory, the following expression must be evaluated,

$$Z_{ij}^{ab} = \sum_{me} t_{im}^{ae} {\cal W}_{mbej}.$$

An efficient approach would first ensure that the summation indices m and e, which are common to both t and ${\cal W}$, are collected together and ordered identically in each. That is, for a given value of each of the ``external'' indices i, j, a, and b, the values of t and ${\cal W}$ for all m and e should be ordered consecutively in memory. If it is possible to store all of t or ${\cal W}$ in core (perhaps along with an additional scratch array), any rearrangment of the data is feasible. It is common, however, for the dimension of such arrays to far exceed the memory capabilities of modern computer hardware, and other schemes must be found.

In 1969, Yoshimine[1,2] developed an efficient, general sorting algorithm for two-electron integrals, but which may be applied to any four-index quantity. In this approach, the values of the array X(pq,rs) to be sorted are assumed to be stored on disk in a random ordering along with the indices p, q, r, and s. As a specific example, consider the ``presorting'' of the atomic-orbital-basis two-electron integrals in the first step of the common integral transformation needed for the construction of many types of correlated wave functions. We require that the integrals (pq|rs) (in standard Mulliken notation) be ordered in an X(pq,rs)``supermatrix'' with the compound indices pq and rs labelling the rows and columns, respectively. An example of the form of such a supermatrix for a calculation involving only four basis functions is illustrated in Fig. 1. Since the two-electron integrals are symmetric with respect to permutation of the bra or ket indices (assuming real basis functions), we wish to store only the unique $p \geq q$ and $r \geq s$ pairs.[3] For example, the compound row index pq is determined from the individual p and qvalues by the equation[4]

$$pq = \left\{ \begin{array}{ll} p (p + 1)/2 + q & \mbox{if $p ... ... q (q + 1)/2 + p & \mbox{if $q \geq p$ } \end{array}\right.$$

If there are N basis functions, the total number of pq or rspairs is N(N+1)/2.[5]

The Yoshimine sort proceeds as follows:

• First, determine how many rows (pq pairs) can be stored in memory at once. Each of these collections of rows is referred to by Yoshimine as a ``core load.'' If we have four basis functions, for example, there will be 4(4+1)/2 = 10 rows/columns in the supermatrix (cf. Fig. 1). Given an unbelievably small machine with only 256 bytes (= 32 eight-byte words) of core memory, we can store three rows in memory simultaneously, and there will be four core loads, the last of which contains only one row of the supermatrix.
• Next, divide the total memory into buffers, each of which corresponds to one core load. In our working example, we would divide the 256 bytes of memory into four buffers of 64 bytes each. Each of these buffers is sufficient to hold only eight integrals at a time.
• Examine each of the integrals stored on disk and place them into the appropriate buffers. We use the indices p and q stored with the integrals to determine to which core load the given integral belongs. For example, if the integral has a p index of 2 and a qindex of 1, it belongs to row 4. In the four basis function example above, this integral would be placed in the second core load.
• The memory buffers are not large enough to hold a complete core load, and once a buffer is filled in memory, its contents are written to disk (in the order in which they arrived) and the memory cleared to allow for more integrals.
• After all the integrals on disk have been distributed to their proper buffers, each core load is read from disk into memory, sorted individually, and written back out to disk. For our working example, the integrals for each core load (which correspond to groups of three consecutive rows of the supermatrix) would be read into memory and sorted ``in place'' such that they have the proper rs pair ordering. The sorted core load may then written back to a final sorted integral file for use in the subsequent integral transformation.

The Yoshimine sorting algorithm is flexible enough to allow for essentially any desired re-ordering of four-index quantities. One need only be able to map individual p, q, r, and s indices to compound row and column indices of the desired supermatrix; otherwise the algorithm remains the same. The resource requirements of this approach are modest: one must be able to store at least one complete row of the supermatrix in memory. For the two-electron integral example above this amounts to only an N² core dependence. The approach can also be extended to arrays with more than four indices with no significant difficulties, though the memory requirements may obviously change. More importantly, however, the Yoshimine sorting algorithm can -- without any adjustment of the working computer code -- take advantage of the large amount of memory available on modern workstations, but still work effectively with the relatively small core sizes found on personal computers.

  

Figure 1: The structure of the X(pq,rs) supermatrix of sorted two-electron integrals for the simple case of four basis functions.