NEW ALGORITHMS OF FAST DIAGONALIZATION OF REAL SYMMETRIC MATRICES.

Author: Yurii F. Sigolaev (St. Petersburg). E-mail: y_sigolaev@hotmail.com

Last revision: July 11, 2007 (The address of page for processor Core 2 Duo: http://www.thesa-store.com/products/)

     The performance of personal computers was dramatically improved during the last decade. Therefore, cardinal revision of the majority of algorithms becomes an urgent problem. In so doing, it is necessary to lift the contradiction between the 32-bit architecture (imposing significant limitations upon the volume of information being processed) and increased processor speed.

      I have developed a new complex of algorithms of diagonalization of dense real symmetric matrices taking into account this factor. This complex (SDIAG) is realized in C and has a number of advantages compared to other known packages realizing diagonalization algorithms:

-          noticeable increase in the computation speed;

-          considerable saving of RAM. The difference (by a factor of 6 when finding all eigenvectors and by a factor of 8 when finding a part of eigenvectors) between all the known high-speed modern diagonalization methods and the algorithms I suggest is so large that my algorithms will apparently remain urgent with the 64-bit architecture also.

      These results were attained thanks to the solution of a basic problem of numerical analysis of tridiagonal matrices, development of new algorithm for vector multiplication and improvement of the block algorithms for packed matrices. Some of the algorithms use the IEEE-arithmetic support and use SSE2 for efficiency. For a comparative analysis, I used Intel MKL package because of the common source of the code, LAPACK package.

      Application of the new algorithm for vector multiplication allows us to note that a pearl in the Intel MKL package, the procedure for the finding full spectrum of symmetric matrix DSYEVD Intel MKL is ineffective (see Sections 1.1, 3.1 and 4).
       I improved the classical Pal-Walker-Kahan's algorithm (section 5).

      The correctness of the SDIAG code was checked within the framework of testing means included in the CLAPACK package (section 8). SDIAG also successfully passed all the tests sent me by Intel.

      The diagonalization program was fairly exhaustively tested with GAMESS quantum-chemical program [12] in which it was included with M.W. Schmidt's permission. The reliability of the results obtained was judged from their coincidence with the results of computation by the PC GAMESS program known for its high speed [13].

      The new results, indicative of a breakthrough in one of the most interesting fields of numerical analysis from both theoretical and applied viewpoints, are also described in [9-11].

Please address all remarks and suggestions to: y_sigolaev@hotmail.com

Tel/Fax: 7(812)2691898

  The program purpose

    A dramatic growth in the performance of personal computers allows a considerable increase in the order of matrices that can be diagonalized within relatively short computation time. However, serious limitations upon the order of matrices being diagonalized are imposed by the 32-bit architecture, because all the modern exact high-speed diagonalization methods require large additional RAM. A primitive way to solve this problem is the use of multiprocessor clusters. I suggest another approach based on adaptation of the existing algorithms to the modern conditions and on the development of new algorithms to overcome the limitations of the 32-bit architecture.

    Note that the algebraic problem of eigenvalues and eigenvectors is a key problem for many scientific and applied tasks.

    The program allows diagonalization of matrices of the following size depending on the RAM size: up to 22000×22000 at 2 GB, (15000×15000 at 1 GB, 11500×11500 at 512 MB, and 7500×7500 at 256 MB, and the diagonalization time is considerably shorter compared to the most known alternative approaches. In my algorithm, the RAM is used only for storing the information only slightly exceeding in size a half of the square matrix. One of the most known and, apparently, the best today is The divide and conquer algorithm allowing diagonalization of matrices of size less than 9000×9000 on a computer with 2 GB RAM. With The divide and conquer algorithm, to find even a part of eigenvectors, the maximal size of matrices being diagonalized decreases to less than 8000×8000 on a computer with 2 GB RAM.

    The algorithms applied to finding eigenvectors of symmetric tridiagonal matrix have strictly quadratic convergence and, at the same time, do not require additional RAM resource.

    The algorithms that I realized for a symmetric tridiagonal matrix develop classical ideas (see Wilkinson's "The Algebraic Eigenvalue Problem"). These algorithms were included for the first time in a quantum-chemical program [1] used for many years for studying the electronic structure of molecules [2-8]. In this study, these algorithms were developed further and tested in the well-known quantum-chemical program GAMESS (see Section 7 and paper [9-11]).

     It should be recognized that the realization of BLAS algorithms in Intel MKL is inefficient (Tables 1, 8). This primarily concerns such functions as DSYMV Intel MKL, DSYR2K Intel MKL, and DGEMM Intel MKL, mainly determining the time required for the tridiagonalization and multiplication of the vectors.

     My program is based in part on the splendid initial code of the LAPACK package which, as follows from the above-said, was subjected to major revision with preservation of all its advantages.

     SDIAG leads to considerable RAM saving in applied tasks (Sections 3.2 and 7; paper [9-11]).

  Comparison results.

     I compared my results with those obtained using the Intel MKL package (Package ID: w_mkl_b_8.0.005 and w_mkl_p_8.1.014).

1.1 My tridiagonalization program operates with packed matrices and therefore requires two times less memory than does DSYTRD Intel MKL (realization of Householder's method for a square matrix with block treatment) and considerably surpasses it in the speed (P4 3.6 GHz processor (Processor 560), Data cache size : L1 16 KB, L2 1024 KB, L3 0 KB, FSB 800 MHz, Motherboard: P5AD2-E Premium, DDR2 533 MHz 2GB, OS Windows XP Professional SP2). The results are compared in Table 1. The difference in the speeds is determined by improved realization of the block methods for packed matrices. (Note: The block size for DSYTRD Intel MKL depends on the size of the apportioned additional RAM. For the matrices under consideration, it is 32).

TABLE 1

1.2 Intel MKL has an option for tridiagonalization of packed matrices, but their speed is extremely low. The comparison is given in Table 2.

TABLE 2

2.1 Algorithms that I realized for a symmetric tridiagonal matrix do not consume additional RAM, have strictly quadratic convergence, and allow separate calculation of eigenvectors of a tridiagonal matrix, irrespective of the distribution of eigenvalues. The calculated eigenvectors are orthogonal within the computer precision.

Consider a matrix W^-_2n+1 for which all the eigenvalues are well separated:

The example was taken from J.H. Wilkinson's monograph "The Algebraic Eigenvalue Problem."

The following results were obtained when finding all eigenvalues and eigenvectors of the W^-_2n+1 matrix using SDIAG:

TABLE 3

2.2 If all the eigenvalues are quasi-degenerated, the calculation time for the total spectrum becomes 2 to 30% longer.

Consider the W⁺_2n+1 matrix:

Almost all the eigenvalues of the W⁺_2n+1 matrix are quasi-degenerated pairwise. Relatively Robust Representations algorithms of finding eigenvalues and eigenvectors of a symmetric tridiagonal matrix, on which much hopes were set, give very bad eigenvectors for the W⁺_2n+1 matrix. The results are compared in Table 4.

TABLE 4

2.3 Consider a symmetric tridiagonal matrix with diagnoal elements equal to 0.1 and nondiagonal elements equal to 2.3e-17. All the eigenvalues of this matrix are quasi-degenerated and equal to 0.1. DSTEGR Intel MKL fails to treat this matrix because of an internal error (info = 2). The results of my calculations coincide in time with the results given in Table 4 to within 5%. The eigenvectors are orthogonal with the computer precision.

3.1  Using of full spectrum of tridiagonal matrix in the completing procedure of Householder method leas to a jump-like decrease in time for calculation of eigenvectors of the parent matrix from the eigenvectors of tridiagonal matrix. This result is achieved owing to application of new algorithm for vector multiplication. A comparison is given in Table 5.

TABLE 5

3.2 Since the eigenvectors of a tridiagonal matrix are calculated separately and independently of the eigenvalue distribution [Section 2.1], they can be found as small blocks, with subsequent multiplication of each block by the matrix obtained by the tridiagonalization and saving of this block on an external carrier in the case of RAM deficiency. This organization of the computations allows significant saving of RAM, so that RAM is used only for storage of the information only slightly exceeding in size a half of the square matrix. Therefore, it becomes feasible to diagonalize matrices of sizes 22000×22000 (2 GB RAM), 15000×15000 (1 GB), 11500×11500 (512 MB), or 7500×7500 (256 MB). The calculation speed is higher than in DORMTR Intel MKL (completing procedure in Householder's method). The results are compared in Table 6. (Note: The block size for DORMTR Intel MKL depends on the size of the apportioned additional RAM. For the matrices under consideration, it is 32).

   In the above-described computation procedure, storage of the eigenvectors is often unnecessary, because in the majority of applied tasks they are an intermediate step in the computation. In this case, the gain in the memory is significant. The RAM is used only for storage of the information only slightly exceeding in size a half of the square matrix (an example of a task of this type is given in Section 7).

TABLE 6

TABLE 7

TABLE 8

4.  Consider a packed symmetric Toplitz matrix:

A_ij = arcsin(1 / (i - j + 1))/(i - j + 1)³ (i , j = 1, 2, ..., n; i >= j)

Application of the algorithms described above in Sections 1.1, 2 and 3.1 to the solution of full algebraic problem for this matrix gave results which shows that the DSYEVD Intel MKL procedure (currently the only a procedure joining calculation speed and accuracy but requiring a huge memory resource) is not effective.

ТАБЛИЦА 10

5.1 Modified procedure for calculating all eigenvalues of a symmetrical tridiagonal matrix, using the QL or QR version of Pal-Walker-Kahan's algorithm with my improvements and allowing determination of the spectrum of the matrix considerably faster than does DSTERF Intel MKL.

    Consider a matrix:

The results are compared in Table 11.

TABLE 11

5.2 Consider a tridiagonal matrix with diagonal elements equal to 2 and nondiagonal elements equal to 1. The eigenvalues of this matrix are equal to 4*sin²[π*k/2*(n + 1)] (k = 1, 2, ..., n). Let us construct vector equal to the difference between the vectors of the exact and approximate solutions. Their norms l_∞ and l₁ characterizing the accuracy of approximate solution show that the SDIAG results have higher precision (Tables 12, 13).

TABLE 12

TABLE 13

6. Below is an example of applying DSTEVD Intel MKL to a matrix from Section 5.2, which demonstrates the cubic nature of The divide and conquer algorithm.

TABLE 14

Note that the oddness of the matrix is insignificant here.

7. Quantum-chemical calculations of large molecules by semiempirical methods (PM3, AM1, etc.) on a single computer, under the conditions of the considerably increased data processing speed in the past decade, are mainly limited by the time and memory resources spent for multiple diagonalization of the Fock matrix. To widen the limits of applicability of a personal computer to quantum-chemical studies of polyatomic compounds, we replaced the diagonalization algorithms included in the original version of the GAMESS program [12] by SDIAG and calculated with this program by the PM3 method a cluster consisting of 14 molecules of fullerene C₆₀ with the parameters of the face-centered crystal lattice of fullerene (840 atoms). The order of the matrices being diagonalized is 3360. Each fullerene molecule has an icosahedral symmerty I_h, and the symmetry of the whole system is decreased to C_i. The symmetry of the system was not used in the calculations. Because of the relatively weak interaction between the fullerene molecules, there are numerous quasi-degenerated energy levels, which makes the testing of the eigenvectors more rigorous. To save RAM, the SDIAG algorithms in the SCF procedure are combined with the improved algorithm of calculating the density matrix, so that the storage of the matrix of eigenvectors of the Fock matrix becomes unnecessary.

   The use of the SDIAG algorithm in this study shortened the time spent for one iteration in the SCF procedure in calculation of the С₈₄₀ cluster to 28.9 s, of which 25.0 s is spent for the diagonalization of the Fock matrix.

    Using in the multiplication procedure the algorithm that multiplies all vectors necessary for construction of density matrix at once reduces the time for one SCF iteration to 27.1 sec of which 23.1 sec is consumed for diagonalization of Fock matrix.

   8. The correctness of the SDIAG code was checked within the framework of testing means included in the CLAPACK package. The corresponding interface has been realized for the SDIAG functions. The results of testing of SDIAG coincide with results of testing of source code of the CLAPACK package:
 

 Tests of the Symmetric Eigenvalue Problem routines
LAPACK VERSION 3.0, released June 30, 1999
The following parameter values will be used:
M: 0 1 2 3 5 20
N: 0 1 2 3 5 20
NB: 1 3 3 3 10
NBMIN: 2 2 2 2 2
NX: 1 0 5 9 1

Relative machine underflow is taken to be .222507-307
Relative machine overflow is taken to be .179769+309
Relative machine precision is taken to be .111022E-15

Routines pass computational tests if test ratio is less than 50.00

DST routines passed the tests of the error exits (147 tests done)

SEP: NB = 1, NBMIN = 2, NX = 1
All tests for DST passed the threshold ( 4662 tests run)
All tests for DST drivers passed the threshold ( 14256 tests run)

SEP: NB = 3, NBMIN = 2, NX = 0
All tests for DST passed the threshold ( 4662 tests run)
All tests for DST drivers passed the threshold ( 14256 tests run)

SEP: NB = 3, NBMIN = 2, NX = 5
All tests for DST passed the threshold ( 4662 tests run)
All tests for DST drivers passed the threshold ( 14256 tests run)

SEP: NB = 3, NBMIN = 2, NX = 9
All tests for DST passed the threshold ( 4662 tests run)
All tests for DST drivers passed the threshold ( 14256 tests run)

SEP: NB = 10, NBMIN = 2, NX = 1
All tests for DST passed the threshold ( 4662 tests run)
All tests for DST drivers passed the threshold ( 14256 tests run)

End of tests
Total time used = 1.53 seconds

References

[1] Yu.F. Sigolaev and S.G. Semenov, Vestn. Leningr. Gos. Univ., 1985, no. 11, pp. 113-114.

[2] Yu.F. Sigolaev, S.G. Semenov, and V.O. Reikhsfel'd, Teor. Exp. Khim., 1983, vol. 19, no. 3, pp. 348-351.

[3] Yu.F. Sigolaev, S.G. Semenov, and V.O. Reikhsfel'd, Teor. Exp. Khim., 1984, vol. 20, no. 4, pp. 482-485.

[4] S.G. Semenov and Yu.F. Sigolaev, Koord. Khim., 1985, vol. 11, no. 12, pp. 1635-1638.

[5] S.G. Semenov, Yu.F. Sigolaev, V.V. Redchenko, and Ya.F. Freimanis, Zh. Obshch. Khim., 1985, vol. 55, no. 2, pp. 401-404.

[6] S.G.Semenov, Yu.F. Sigolaev. and S. M. Shevchenko, Croat. Chem. Acta, 1988, vol. 61, no. 1, pp. 73-80.

[7] S.G. Semenov and Yu.F. Sigolaev, Koord. Khim., 1988, vol. 14, no.12, pp. 1636-1640.

[8] Yu.F. Sigolaev, S.G. Semenov, V.O. Reikhsfel'd, and V.K. Bel'skii, Teor. Exp. Khim.,1990, vol. 26, no. 2, pp. 221-224.

[9] S.G. Semenov and Yu.F. Sigolaev, Zh. Obshch. Khim., 2006, vol. 76, no. 6, pp. 1006-1009.

[10] S.G. Semenov and Yu.F. Sigolaev, Zh. Obshch. Khim., 2006, vol. 76, no. 11, pp. 1732-1737.

[11] S.G. Semenov and Yu.F. Sigolaev, Zh. Obshch. Khim., 2007, vol. 77, no. 5, pp. 727-732.

[12] M.W. Schmidt, K.K. Baldridge, J.A. Boatz et al., J. Comput. Chem., 1993, vol. 14, no. 9, p. 1347.

[13] A.A. Granovsky. http: // classic.chem.msu.su/gran/gamess/index.html.