_b_e_g_i_n _c_o_m_m_e_n_t eigenvalues of a real symmetric matrix
by the QR method. Algorithm 253, P.A. Businger, CACM 8 (1965) 217;
_i_n_t_e_g_e_r n;
n:= read;
_b_e_g_i_n _i_n_t_e_g_e_r i,j;
_r_e_a_l _a_r_r_a_y a[1:n,1:n];
_p_r_o_c_e_d_u_r_e symmetric QR1(n,g);
_v_a_l_u_e n; _i_n_t_e_g_e_r n; _a_r_r_a_y g;
_c_o_m_m_e_n_t 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;
_b_e_g_i_n
_r_e_a_l _p_r_o_c_e_d_u_r_e sum(i,m,n,a);
_v_a_l_u_e m,n; _i_n_t_e_g_e_r i,m,n; _r_e_a_l a;
_b_e_g_i_n _r_e_a_l s; s:= 0;
_f_o_r i:= m _s_t_e_p 1 _u_n_t_i_l n _d_o s:= s + a;
sum:= s
_e_n_d sum;
_r_e_a_l _p_r_o_c_e_d_u_r_e max (a,b);
_v_a_l_u_e a,b; _r_e_a_l a,b;
max:= _i_f a > b _t_h_e_n a _e_l_s_e b;
_p_r_o_c_e_d_u_r_e Housholder tridiagonalization 1(n,g,a,bq,norm);
_v_a_l_u_e n; _i_n_t_e_g_e_r n; _a_r_r_a_y g,a,bq; _r_e_a_l norm;
_c_o_m_m_e_n_t nonlocal real procedure sum, max;
_c_o_m_m_e_n_t 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;
_b_e_g_i_n _i_n_t_e_g_e_r i,j,k;
_r_e_a_l t,absb,alpha,beta,gamma,sigma;
_a_r_r_a_y p[2:n];
norm:= absb:= 0;
_f_o_r k:= 1 _s_t_e_p 1 _u_n_t_i_l n - 2 _d_o
_b_e_g_i_n a[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);
_i_f sigma |= 0 _t_h_e_n
_b_e_g_i_n alpha:= g[k+1,k];
beta:= _i_f alpha < 0 _t_h_e_n absb _e_l_s_e - absb;
gamma:= 1 / (sigma-alpha*beta); g[k+1,k]:= alpha - beta;
_f_o_r i:= k + 1 _s_t_e_p 1 _u_n_t_i_l n _d_o
p[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]);
_f_o_r i:= k + 1 _s_t_e_p 1 _u_n_t_i_l n _d_o
p[i]:= p[i] - t*g[i,k];
_f_o_r i:= k + 1 _s_t_e_p 1 _u_n_t_i_l n _d_o
_f_o_r j:= k + 1 _s_t_e_p 1 _u_n_t_i_l i _d_o
g[i,j]:= g[i,j] - g[i,k]*p[j] - p[i]*g[j,k]
_e_n_d
_e_n_d k;
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]))
_e_n_d Housholder tridiagonalization 1;
_i_n_t_e_g_e_r i,k,m,m1;
_r_e_a_l norm,epsq,lambda,mu,sq1,sq2,u,pq,gamma,t;
_a_r_r_a_y a[1:n],bq[0:n-1];
Housholder tridiagonalization 1(n,g,a,bq,norm);
epsq:= 2.25%-22*norm|^2;
_c_o_m_m_e_n_t 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: _i_f m = 0
_t_h_e_n _g_o_t_o return _e_l_s_e i:= k:= m1:= m - 1;
bq[0]:= 0;
_i_f bq[k] _< epsq _t_h_e_n
_b_e_g_i_n g[m,m]:= a[m]; mu:= 0; m:= k;
_g_o_t_o inspect
_e_n_d;
_f_o_r i:= i - 1 _w_h_i_l_e bq[i] > epsq _d_o k:= i;
_i_f k = m1 _t_h_e_n
_b_e_g_i_n _c_o_m_m_e_n_t 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*(_i_f sq1_>0 _t_h_e_n sq1+sq2 _e_l_s_e sq1-sq2);
g[m1,m1]:= lambda; g[m,m]:= mu / lambda;
mu:= 0; m:= m - 2; _g_o_t_o inspect
_e_n_d;
lambda:= _i_f abs(a[m]-mu) < 0.5*abs(a[m])
_t_h_e_n a[m] + 0.5*sqrt(bq[m1]) _e_l_s_e 0.0;
mu:= a[m]; sq1:= sq2:= u:= 0;
_f_o_r i:= k _s_t_e_p 1 _u_n_t_i_l m1 _d_o
_b_e_g_i_n _c_o_m_m_e_n_t shortcut single QR iteration;
gamma:= a[i] - lambda - u;
pq:= _i_f sq1 |= 1
_t_h_e_n gamma|^2/(1-sq1) _e_l_s_e (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 + lambda
_e_n_d i;
gamma:= a[m] - lambda - u;
bq[m1]:= sq1 *
(_i_f sq1|=1 _t_h_e_n gamma|^2/(1-sq1) _e_l_s_e (1-sq2)*bq[m1]);
a[m]:= gamma + lambda; _g_o_t_o inspect;
return:
_e_n_d symmetric QR 1;
_f_o_r i:= 1 _s_t_e_p 1 _u_n_t_i_l n _d_o
_f_o_r j:= 1 _s_t_e_p 1 _u_n_t_i_l i _d_o
a[i,j]:= read;
symmetric QR 1(n,a);
_f_o_r i:= 1 _s_t_e_p 1 _u_n_t_i_l n _d_o
_b_e_g_i_n NLCR; ABSFIXT(2,0,i); PRINT(a[i,i]) _e_n_d
_e_n_d
_e_n_d