// ************************************************************************** // // // // eses eses // // eses eses // // eses eseses esesese eses Embedded Systems Group // // ese ese ese ese ese // // ese eseseses eseseses ese Department of Computer Science // // eses eses ese eses // // eses eseses eseseses eses University of Kaiserslautern // // eses eses // // // // ************************************************************************** // // The algorithm below improves module EigenvaluesQR in that shift operations // // are used to speed-up the QR-iteration. The idea is as follows: instead of // // using A_0 := A and A_{k+1} := R_k*Q_k where A_k = Q_k * R_k is obtained by // // the QR-decomposition of A_k, we use A_0 := A and A_{k+1} := R_k*Q_k+p_k*E // // where A_k-p_k*E=Q_k*R_k is obtained by the QR-decomposition of A_k-p_k*E. // // It is then easily seen that A_{k+1} = Q_k^T*A_k*Q_k so that all A_k have // // the same eigenvalues, regardless which value is chosen for p_k. It can be // // proved that the convergence of A[i][j]_k is linear with |e[i]/e[j]|^k. Now,// // note that A_k-p_k*E has the eigenvalues e[i]-p_k so that A[i][j]_k will // // converge as |(e[i]-p_k)/(e[j]-p_k)|^k. Thus, one should choose p_k close to// // e[i+1], but since the latter is not yet known, one has to guess. // // In particular, one often chooses A[i][i], but it is known that this does // // not work with complex roots. Wilkinson therefore proposed to use the // // eigenvalue of the 2x2 matrix [[A[i][i],A[i][i+1]],[A[i][i+1],A[i+1][i+1]]] // // that is closer to A[i+1][i+1] (note that this eigenvalue may be complex). // // ************************************************************************** // macro N = 5; // numeric precision macro rabs(x) = (x<0 ? -x : +x); macro sqrt(x) = exp(x,0.5); macro eps = 1.0e-10; module EigenvaluesShiftQR([N][N]real ?m,a,[N]real e,event !rdy) { [N][N]real q,r; // matrices q,r of QR-decomposition real p; // shift parameter nat{N+1} n; // current size of the matrix (due to deflation) event QRrdy; // whether QR-decomposition terminated // initialize n = N; for(i=0..N-1) for(j=0..N-1) a[i][j] = m[i][j]; do { do { // determine shift p: note that this simple approximation to e[N-1] // works well in case the latter is a real number; otherwise the // QR-iteration will not converge p = a[n-1][n-1]; // shift matrix a for(i=0..N-1) next(a[i][i]) = a[i][i]-p; pause; // decompose shifted matrix A_k-p*E into factors Q_k and R_k DecompQR(a,q,r,QRrdy); // compute A_{k+1} := R_k*Q_k+p_k*E and check for convergence for(i=0..N-1) { for(j=0..N-1) { let(p0 = (i==j ? p : 0.0)) next(a[i][j]) = p0 + sum(k=0..N-1) (r[i][k]*q[k][j]); } } pause; } while(rabs(a[n-1][n-2])>eps); next(n) = n-1; // eigenvalue e[n-1] has been found, we can deflate the matrix e[n-1] = a[n-1][n-1]; for(i=0..N-1) { next(a[i][n-1]) = 0.0; next(a[n-1][i]) = 0.0; } pause; } while(n>1); e[n-1] = a[n-1][n-1]; emit(rdy); } drivenby m0 { m[0] = [1.0,2.0,3.0,4.0,5.0]; m[1] = [2.0,3.0,4.0,5.0,1.0]; m[2] = [3.0,4.0,5.0,1.0,2.0]; m[3] = [4.0,5.0,1.0,2.0,3.0]; m[4] = [5.0,1.0,2.0,3.0,4.0]; await(rdy); } drivenby m1 { // note: x^5 - 15x^4 + 85x^3 - 225x^2 + 274x - 120 = (x-1)(x-2)(x-3)(x-4)(x-5) m[0] = [ 15.0, -85.0, 225.0, -274.0, 120.0]; m[1] = [ 1.0, 0.0, 0.0, 0.0, 0.0]; m[2] = [ 0.0, 1.0, 0.0, 0.0, 0.0]; m[3] = [ 0.0, 0.0, 1.0, 0.0, 0.0]; m[4] = [ 0.0, 0.0, 0.0, 1.0, 0.0]; await(rdy); } drivenby m2 { // note: x^5 - 2x^4 - 3x^3 + 6x^2 + 2x - 4 = (x-2)(x-1)(x+1)(x-sqrt(2))(x+sqrt(2)) m[0] = [ 2.0, 3.0, -6.0, -2.0, 4.0]; m[1] = [ 1.0, 0.0, 0.0, 0.0, 0.0]; m[2] = [ 0.0, 1.0, 0.0, 0.0, 0.0]; m[3] = [ 0.0, 0.0, 1.0, 0.0, 0.0]; m[4] = [ 0.0, 0.0, 0.0, 1.0, 0.0]; await(rdy); } // m3 does not converge since there are complex eigenvalues //drivenby m3 { // note: x^5+6x^4+14x^3+24x^2+33x+18 = (x+3)(x+2)(x+1)(x^2+3) // m[0] = [ -6.0, -14.0, -24.0, -33.0, -18.0]; // m[1] = [ 1.0, 0.0, 0.0, 0.0, 0.0]; // m[2] = [ 0.0, 1.0, 0.0, 0.0, 0.0]; // m[3] = [ 0.0, 0.0, 1.0, 0.0, 0.0]; // m[4] = [ 0.0, 0.0, 0.0, 1.0, 0.0]; // await(rdy); //} // m4 does not converge since there are complex eigenvalues //drivenby m4 { // note: x^5 - 2x^4 + 2x^3 - 4x^2 - 3x + 6 = (x-2)(x-1)(x+1)(x^2+3) // m[0] = [ 2.0, -2.0, 4.0, 3.0, -6.0]; // m[1] = [ 1.0, 0.0, 0.0, 0.0, 0.0]; // m[2] = [ 0.0, 1.0, 0.0, 0.0, 0.0]; // m[3] = [ 0.0, 0.0, 1.0, 0.0, 0.0]; // m[4] = [ 0.0, 0.0, 0.0, 1.0, 0.0]; // await(rdy); //}