(HP71B) Eigenvalues of Symmetric Real Matrix

[Albert Chan]

Algorithm from HP15C Advanced Function Handbook, page 125
This version returned both eigenvectors and eigenvalues of symmetric real matrix

10 DESTROY ALL @ OPTION BASE 1 @ INPUT "N=";N
20 DIM A(N,N),B(N,N),L(N,N),Q(N,N),I0(N,N) @ MAT I0=IDN @ MAT Q=I0
30 MAT INPUT A @ MAT L=TRN(A) @ MAT A=A+L @ MAT A=(.5)*A @ MAT L=A
40 DEF FNF(N,D)=N/(D+(1-2*(D<0))*SQRT(N*N+D*D)) ! = tan(atan(N/D)/2)
100 MAT B=ZER
110 FOR I=2 to N @ FOR J=1 to I-1
120 X=L(I,J) @ IF X THEN B(I,J)=FNF(FNF(2*X,L(I,I)-L(J,J)),1)
130 NEXT J @ NEXT I
140 MAT L=TRN(B) @ MAT B=B-L @ DISP "FNORM(B)=";FNORM(B)
150 MAT L=B+I0 @ MAT L=INV(L) @ MAT L=L+L @ MAT L=L-I0 @ MAT Q=Q*L
160 DISP "Q=" @ MAT DISP Q
170 MAT L=TRN(Q) @ MAT L=L*A @ MAT L=L*Q
180 DISP "L=" @ MAT DISP L

>run
N=3
A(1,1)? 0,1,2
A(2,1)? 1,2,3
A(3,1)? 2,3,4

Code:

FNORM(B)= .605552349352
Q=
 .8662568603          .227657549016        .444714619016
-.444714619012        .75699314672         .478738220168
-.227657549018       -.612481359872        .75699314673
L=
-.349003969447       -.425695598709       -1.12609254611
-.425695598709       -.348327850394       -4.47666317212E-2
-1.12609254611       -4.47666317247E-2     6.69733181992
>run 100
FNORM(B)= .595643959403
Q=
 .824769255982       -.469983482558        .314437912022
 .285388154078        .826016169657        .486056466843
-.488169310695       -.311147575236        .815400460218
L=
-.859525935837        9.56986605065E-2     .200156808749
 9.56986605065E-2     1.82880814699E-2     .42420712414
 .200156808747        .424207124136        6.84123785436
>run 100
FNORM(B)= 8.95738317424E-2
Q=
 .86332901543        -.398599604012        .309485648795
 .183632804446        .819374127488        .543051592685
-.470044683262       -.412000479214        .78058542189
L=
-.872825507534       -9.30648141545E-3    -2.14243725402E-2
-9.30648141545E-3    -8.11920389207E-5     1.09063583269E-2
-2.14243725373E-2     .010906358327        6.87290669953
>run 100
FNORM(B)= 7.86910670518E-3
Q=
 .859879422537       -.408274692867        .306462320485
 .193856724533        .816499388716        .543827471319
-.472257291179       -.408216270356        .78123781753
L=
-.872983344667        3.2546038359E-5      5.1700251917E-5
 3.2546038358E-5      1.1264547796E-9      1.26811006234E-4
 5.170025078E-5       1.2681100214E-4      6.87298334337
>run 100
FNORM(B)= 2.97899107959E-5
Q=
 .859892597624       -.408248289526        .30646052708
 .193822654531        .816496581206        .54384383001
-.472247286024       -.408248290825        .781227133338
L=
-.872983346221       -8.99996E-10         -2.063351E-9
-8.99996266439E-10   -7.5924E-19           1.07648489092E-9
-2.05974E-9           1.07239E-9           6.87298334594
>run 100
FNORM(B)= 7.60886994417E-10
Q=
 .859892597285       -.40824829046         .306460526788
 .193822655517        .816496580921        .543843830086
-.472247286237       -.40824829046         .7812271334
L=
-.872983346221       -1.69E-13            -3.663E-12
-1.69203950392E-13    2.E-24               3.73782958716E-12
 1.66E-12             3.73E-12             6.87298334595

Q^T A Q = L      → A = Q L Q^T

Werner recently PM me on how this work, but I have no clue.
Anyone knows how it work, or reference (besides AFH) to where algorithm comes from?