begincomment eigenvalues of a real symmetric matrix
by the QR method. Algorithm 253, P.A. Businger, CACM 8 (1965) 217;integerN;N:=READ;beginintegerI,J;realarrayA[1:N,1:N];procedureSYMMETRICQR1(N,G);valueN;integerN;arrayG;comment uses Housholders.s method and the QR algorithm to
find all n eigenvalues of the real symmetric matrix whose lower
triangular part is given in the array g[1:n,1:n]. The computed
eigenvalues are stored as the diagonal elements g[i,i]. The
original contents of the lower triangular part of g are lost during
the computation whereas the strictly upper triagular part of g
is left untouched;beginrealprocedureSUM(I,M,N,A);valueM,N;integerI,M,N;realA;beginrealS;S:=0;forI:=Mstep1untilNdoS:=S+A;SUM:=SendSUM;realprocedureMAX(A,B);valueA,B;realA,B;MAX:=ifA>BthenAelseB;procedureHOUSHOLDERTRIDIAGONALIZATION1(N,G,A,BQ,NORM);valueN;integerN;arrayG,A,BQ;realNORM;comment nonlocal real procedure sum, max;comment reduces the given real symmetric n by n matrix g
to tridiagonal form using n - 2 elementary orthogonal trans-
formations (I-2ww.) = (I-gamma uu.). Only the lower tri
angular part of g need be given. The diagonal elements and
the squares of the subdiagonal elements of the reduced matrix
are stored in a[1:n] and bq[1:n-1] respectively. norm is set
equal to the infinity norm of the reduced matrix. The columns
of the strictly lower triagular part of g are replaced by the
nonzero portions of the vectors u;beginintegerI,J,K;realT,ABSB,ALPHA,BETA,GAMMA,SIGMA;arrayP[2:N];NORM:=ABSB:=0;forK:=1step1untilN-2dobeginA[K]:=G[K,K];SIGMA:=BQ[K]:=SUM(I,K+1,N,G[I,K]|^2);T:=ABSB+ABS(A[K]);ABSB:=SQRT(SIGMA);NORM:=MAX(NORM,T+ABSB);ifSIGMA|=0thenbeginALPHA:=G[K+1,K];BETA:=ifALPHA<0thenABSBelse-ABSB;GAMMA:=1/(SIGMA-ALPHA*BETA);G[K+1,K]:=ALPHA-BETA;forI:=K+1step1untilNdoP[I]:=GAMMA*(SUM(J,K+1,I,G[I,J]*G[J,K])+SUM(J,I+1,N,G[J,I]*G[J,K]));T:=0.5*GAMMA*SUM(I,K+1,N,G[I,K]*P[I]);forI:=K+1step1untilNdoP[I]:=P[I]-T*G[I,K];forI:=K+1step1untilNdoforJ:=K+1step1untilIdoG[I,J]:=G[I,J]-G[I,K]*P[J]-P[I]*G[J,K]endendK;A[N-1]:=G[N-1,N-1];BQ[N-1]:=G[N,N-1]|^2;A[N]:=G[N,N];T:=ABS(G[N,N-1]);NORM:=MAX(NORM,ABSB+ABS(A[N-1])+T);NORM:=MAX(NORM,T+ABS(A[N]))endHOUSHOLDERTRIDIAGONALIZATION1;integerI,K,M,M1;realNORM,EPSQ,LAMBDA,MU,SQ1,SQ2,U,PQ,GAMMA,T;arrayA[1:N],BQ[0:N-1];HOUSHOLDERTRIDIAGONALIZATION1(N,G,A,BQ,NORM);EPSQ:=2.25%-22*NORM|^2;comment The tollerance used in the QR iteration depends
on the square of the relative machine precision. Here 2.25%-22
is used which is appropriate for a machine with a 36-bit
mantissa;MU:=0;M:=N;INSPECT:ifM=0thengotoRETURNelseI:=K:=M1:=M-1;BQ[0]:=0;ifBQ[K]<EPSQthenbeginG[M,M]:=A[M];MU:=0;M:=K;gotoINSPECTend;forI:=I-1whileBQ[I]>EPSQdoK:=I;ifK=M1thenbegincomment treat 2 * 2 block separately;MU:=A[M1]*A[M]-BQ[M1];SQ1:=A[M1]+A[M];SQ2:=SQRT((A[M1]-A[M])|^2+4*BQ[M1]);LAMBDA:=0.5*(ifSQ1>0thenSQ1+SQ2elseSQ1-SQ2);G[M1,M1]:=LAMBDA;G[M,M]:=MU/LAMBDA;MU:=0;M:=M-2;gotoINSPECTend;LAMBDA:=ifABS(A[M]-MU)<0.5*ABS(A[M])thenA[M]+0.5*SQRT(BQ[M1])else0.0;MU:=A[M];SQ1:=SQ2:=U:=0;forI:=Kstep1untilM1dobegincomment shortcut single QR iteration;GAMMA:=A[I]-LAMBDA-U;PQ:=ifSQ1|=1thenGAMMA|^2/(1-SQ1)else(1-SQ2)*BQ[I-1];T:=PQ+BQ[I];BQ[I-1]:=SQ1*T;SQ2:=SQ1;SQ1:=BQ[I]/T;U:=SQ1*(GAMMA+A[I+1]-LAMBDA);A[I]:=GAMMA+U+LAMBDAendI;GAMMA:=A[M]-LAMBDA-U;BQ[M1]:=SQ1*(ifSQ1|=1thenGAMMA|^2/(1-SQ1)else(1-SQ2)*BQ[M1]);A[M]:=GAMMA+LAMBDA;gotoINSPECT;RETURN:endSYMMETRICQR1;forI:=1step1untilNdoforJ:=1step1untilIdoA[I,J]:=READ;SYMMETRICQR1(N,A);forI:=1step1untilNdobeginNLCR;ABSFIXT(2,0,I);PRINT(A[I,I])endendend