// ************************************************************************** // // // // 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 approximates the eigenvalues of a given NxN matrix by // // means of the QR-decomposition: Given a matrix A, we define the iteration // // 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. Since Q_k is orthogonal, its inverse matrix is its// // transposed matrix, and therefore A_{k+1} = (Q_k^t*Q_k)*(R_k*Q_k) holds, // // which is equivalent to A_{k+1} = Q_k^t*A_k*Q_k. For this reason, all A_k // // have the same eigenvalues, since e*x = A_{k+1}*x = Q_k^t*A_k*Q_k*x, i.e. // // equivalent to e*(Q_k*x) = A_k*(Q_k*x). In case, the sequence A_k converges // // to a matrix A', it then follows that also Q_k and R_k converge to matrices // // Q' and R' so that A' = Q' * R' = R' * Q' holds, where Q' is orthogonal and // // R' is an upper triangular matrix. A' is the a Schur matrix, i.e. a matrix // // whose entries are all zero with the exception of diagonal elements A[i][i] // // and some of the elements A[i][i+1] and A[i+1][i] so that the diagonal of A // // consists of the eigenvalues and 2x2 matrices that represent two complex // // eigenvalues. The simple QR-iteration as implemented here does however not // // converge if there are complex eigenvalues, and if there are eigenvalues of // // the same absolute value. To solve the latter, we would have to deflate the // // matrix, and to deal with complex eigenvalues, we have to use Wilkinson // // shifts with complex numbers (or use the double step QR-iteration). // // ************************************************************************** // macro N = 5; // numeric precision macro rabs(x) = (x<0 ? -x : +x); macro eps = 1.0e-10; module EigenvaluesQR([N][N]real ?m,a,event !rdy) { [N][N]real q,r; event QRrdy; event iterate; // initialize for(i=0..N-1) for(j=0..N-1) a[i][j] = m[i][j]; do { // decompose A_k = Q_k * R_k DecompQR(a,q,r,QRrdy); // determine A_{k+1} = R_k * Q_k and check convergence for(i=0..N-1) { for(j=0..N-1) { let(new_aij = sum(k=0..N-1) (r[i][k]*q[k][j])) { next(a[i][j]) = new_aij; if(i==j & rabs(a[i][j]-new_aij)>eps) next(iterate) = true; } } } pause; next(iterate) = false; } while(iterate); emit(rdy); } drivenby { 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 { // 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); }