// ************************************************************************** // // // // 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 EigenvaluesShiftQR in that two shift // // operations are merged (to obtain the so-called double shift QR method). // // This is very important to achieve convergence for matrices with complex // // eigenvalues, since the shift p = a[N-1][N-1] used in EigenvaluesShiftQR // // will keep all entries of the matrices Q_k,R_k and A_k in the real numbers. // // Thus, it cannot converge for matrices with complex eigenvalues. One way to // // deal with this problem is to work with complex numbers using Wilkinson's // // shift as shown in module EigenvaluesQRC, another way is to perform a double// // shift with the eigenvalues of the 2x2 matrix [[A[N-2][N-2],A[N-2][N-1]], // // [A[N-2][N-1],A[N-1][N-1]]]. It is interesting that these computation can // // be done with real numbers: Starting with matrix A_k, we consider 2x2 matrix// // in the lower right corner and compute its eigenvalues e1 and e2 (that may // // be complex numbers). If these were used as shifts in QR-steps, we would // // compute B_k := R_k*Q_k + e1*E where Q_k*R_k := A_k - e1*E; and in the next // // step, A_k+1 := R_k'*Q_k' + e2*E where Q_k'*R_k' := B_k - e2*E. Thus, we // // have Q_k'*R_k' = R_k*Q_k + (e1-e2)*E, and by multiplying with Q_k left and // // R_k right, we have Q_k*Q_k'*R_k'*R_k = Q_k*R_k*Q_k*R_k + (e1-e2)*Q_k*R_k, // // i.e., Q_k*Q_k'*R_k'*R_k = Q_k*R_k*(Q_k*R_k + (e1-e2)*E). By definition of // // Q_k and R_k, we have Q_k*Q_k'*R_k'*R_k = (A_k - e1*E) * (A_k - e2*E), so // // that (Q_k*Q_k')*(R_k'*R_k) can also be obtained as the QR-decomposition of // // (A_k-e1*E)*(A_k-e2*E) = A_k*A_k - 2*Re(e1)*A_k + |e_1|^2*E, which is a // // real valued matrix. Considering the 2x2 matrix [[A[N-2][N-2],A[N-2][N-1]], // // [A[N-2][N-1],A[N-1][N-1]]] =:G, we can see that the matrix can be computed // // as follows: (A_k-e1*E)*(A_k-e2*E) = A_k*A_k - trace(G)*A_k + det(G)*E. // // ************************************************************************** // 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 EigenvaluesShift2QR([N][N]real ?m,a,[N]real e,event !rdy) { [N][N]real b,q,r; // matrices q,r of QR-decomposition real trcG,detG; // to determine double shift 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 { // compute trace and determinant of corner matrix for double shift trcG = a[n-2][n-2] + a[n-1][n-1]; detG = a[n-2][n-2] * a[n-1][n-1] - a[n-2][n-1] * a[n-1][n-2]; // compute A_k*A_k - trace(G)*A_k + det(G)*E in matrix b for(i=0..N-1) for(j=0..N-1) next(b[i][j]) = detG - trcG*a[i][j] + sum(k=0..N-1) (a[i][k]*a[k][j]); pause; // decompose double shifted matrix b into factors Q_k and R_k DecompQR(b,q,r,QRrdy); // compute A_{k+1} := Q_k^T * A_k * Q_k and check for convergence for(i=0..N-1) { for(j=0..N-1) { b[i][j] = sum(k=0..N-1) (a[i][k]*q[k][j]); next(a[i][j]) = sum(k=0..N-1) (q[k][i]*b[k][j]); } } pause; // the following termination condition must be adapted for complex roots } 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); } //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); //} //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); //}