Jacobi matrix diagonalization problems

[Spathi]

I need to implement a routine for finding the eigenvalues of a symmetrical matrix (for computing the vibrational frequencies from Hessian). I had already implemented a Jacobi diagonalization algorithm, and in most cases it works properly, but sometimes it crashes. In particular, it crashes when trying to diagonalize this matrix:

Spoiler

8.1285326168826E-6 2.38306206017885E-10 3.47264450495782E-9 -1.86224307187575E-6 2.07291390876557E-7 7.10382787255987E-7 -2.11716580785851E-6 -3.56347908891966E-7 1.55354671458546E-6 -2.48362013270562E-6 1.48849726541923E-7 -2.26671204429494E-6
1.64647438695376E-10 8.12125705093878E-6 -4.83879107332822E-9 2.07007814274288E-7 -2.43912835431489E-6 -2.20280352385421E-6 -3.5557217456388E-7 -2.18989978010799E-6 1.71867197602653E-6 1.48415293298861E-7 -1.82821569427162E-6 4.88002290386578E-7
3.08290296016882E-9 -5.67968728244711E-9 1.57935511836939E-5 4.36226239778081E-7 -1.35265046707309E-6 -4.18603775487232E-6 9.5489784625453E-7 1.05668195539855E-6 -4.18499559959507E-6 -1.39359095947773E-6 3.00508536567709E-7 -4.186569525721E-6
-1.86325688032289E-6 2.04241612048088E-7 4.34344712579062E-7 1.95418631395656E-6 -8.57916701225296E-7 -3.32864408651317E-7 3.65683624492258E-7 -4.96914771976087E-7 -4.2739549957064E-7 -8.38376844511562E-7 1.1924461623516E-6 4.14914896208524E-7
2.0477554872122E-7 -2.4392716938814E-6 -1.34967475345946E-6 -8.58307453234214E-7 4.34465280449163E-6 1.03477945772584E-6 -4.80608873011829E-7 -1.80435674114853E-6 -1.1666024229946E-7 1.17609783444138E-6 -6.00820618891346E-7 1.5501292986552E-7
7.04334309042532E-7 -2.18678919034304E-6 -4.18032436129733E-6 -3.32246982707501E-7 1.03281029104992E-6 1.76480075851612E-6 -3.42186860981594E-7 2.79257875134355E-7 7.79465812236145E-7 1.14414941996306E-7 4.26659960913203E-7 7.79549641714165E-7
-2.11341141715773E-6 -3.49304912451842E-7 9.51800904867838E-7 3.6559452924397E-7 -4.81095952517078E-7 -3.41987890148466E-7 3.00059039784564E-6 1.46255579597976E-6 -7.28812546407063E-7 -1.68580035591889E-6 -7.03732310446134E-7 3.1402312587757E-7
-3.49855212953954E-7 -2.18901866188688E-6 1.05737118658894E-6 -4.97364326761096E-7 -1.80496787787468E-6 2.80150182140227E-7 1.46280118580038E-6 3.29862132839406E-6 -8.08575833376958E-7 -6.87261145055025E-7 2.46834769007884E-7 -3.12305420377653E-7
1.54197784334517E-6 1.70969638443062E-6 -4.18002782799627E-6 -4.27193051892312E-7 -1.16247434667037E-7 7.79477011829923E-7 -7.27167651347098E-7 -8.06458161925154E-7 1.7645343077223E-6 -7.16748790017276E-8 -4.36694156402397E-7 7.79568534026464E-7
-2.47942150730988E-6 1.46386123432345E-7 -1.39233422203437E-6 -8.3858606749021E-7 1.17620571629039E-6 1.13903039990619E-7 -1.68624276791078E-6 -6.87052324017123E-7 -7.23378475087464E-8 4.49621413992641E-6 -6.05533637041756E-7 1.06549478105954E-6
1.46444171529516E-7 -1.82595118756624E-6 2.98559889917643E-7 1.19269426645613E-6 -6.00847959409068E-7 4.27271934582779E-7 -7.03647423852404E-7 2.47020006680735E-7 -4.36998861762308E-7 -6.05474561252034E-7 1.80565742960205E-6 -2.27667180665179E-7
-2.25170804837726E-6 4.8313892479228E-7 -4.18240539470707E-6 4.14392732184937E-7 1.55446860561168E-7 7.79638972611766E-7 3.13251631859053E-7 -3.12517837456977E-7 7.79482822703453E-7 1.06271288410322E-6 -2.27083836082802E-7 1.76634810824549E-6


(12 values in each line, i.e. each row in the matrix consists of 3 lines here)

I implemented this algorithm quite long ago and partly forgot how it worked. At each iteration, it finds the biggest non-diagonal element of the matrix, computes some arctangent from the values by its position; then it creates a rotational matrix – firstly it is an identity matrix, but then it replaces 4 digits in it (2 at the diagonal) with sinuses and cosines of this arctangent. Then the source matrix is multiplied by this rotational matrix (more exactly, two “symmetrical” multiplications are made). The iterations are repeated until the non-diagonal elements are sufficiently small.
I can write more details about this algorithm, if this is necessary; it is interesting for me to understand the principle, how all these things work.

 

[marcusl]

Use Matlab, a routine from Lapack, or any number of tried and true packages. There’s no sense in reinventing the wheel. Also check the condition number of your matrix (use the cond command in Matlab). Your matrix may be ill-conditioned.

 

[Spathi]

Use Matlab, a routine from Lapack, or any number of tried and true packages. There’s no sense in reinventing the wheel. Also check the condition number of your matrix (use the cond command in Matlab). Your matrix may be ill-conditioned.

I am creating software with Delphi, so I have to reinvent the wheels.

 

[renormalize]

Your matrix may be ill-conditioned.

Importing the matrix into Mathematica gives the ratio of the largest-to-smallest singular values to be nearly 10^8. @Spathi, is your algorithm designed to handle matrices that are this ill-conditioned?

 

[Spathi]

@Spathi, is your algorithm designed to handle matrices that are this ill-conditioned?

I don't know. Evidently the matrix in the op is correctly diagonalized by Gaussian software (I extracted this matrix from Gaussian .fch file with vibrational frequencies computation).

 

[Spathi]

I want to understand the principles of this calculations. Is it correct that for example for the 12*12 matrix in the op, we can say that initially we had 12 vectors in 12-dimensional space, and then we find a 12*12 rotation matrix that rotates these vectors in a certain way?
Maybe I will be able understand this using the example of calculating the moments of inertia of a molecule:

Here, firstly a symmetric 3*3 matrix is initially calculated with elements like

$$ \sum_{n} m_n y_n^2 z_n^2 $$


Then this matrix is diagonalized; i.e., maybe this matrix is simply rotated by multiplying by a 3*3 rotation matrix? And the rotation matrix is determined by the x,y,z vectors at the image.