// ************************************************************************** // // // // 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 is a version of EigenvaluesShiftQR that// // works on complex numbers. // // ************************************************************************** // macro N = 5; // -------------------------------------------------------------------------- // // numeric precision // -------------------------------------------------------------------------- // macro eps = 1.0e-8; // -------------------------------------------------------------------------- // // macros on real numbers // -------------------------------------------------------------------------- // macro rsign(x) = (x<0.0 ? -1.0 : +1.0); macro rabs(x) = (x<0.0 ? -x : +x); macro rroot(x) = exp(x,0.5); // -------------------------------------------------------------------------- // // arithmetic on complex numbers // -------------------------------------------------------------------------- // macro cadd(x,y) = (x.0+y.0,x.1+y.1); macro csub(x,y) = (x.0-y.0,x.1-y.1); macro cmul(x,y) = (x.0*y.0-x.1*y.1,x.1*y.0+x.0*y.1); macro cdiv(x,y) = ((x.0*y.0+x.1*y.1)/(y.0*y.0+y.1*y.1), (x.1*y.0-x.0*y.1)/(y.0*y.0+y.1*y.1)); macro sqabs(x) = (x.1*x.1 + x.0*x.0); macro cabs(x) = rroot(x.1*x.1 + x.0*x.0); macro croot(x) = (rsign(x.1) * rroot(0.5 * (cabs(x) + x.0)), rroot(0.5 * (cabs(x) - x.0))); macro Re(x) = x.0; macro Im(x) = x.1; macro unit = (1.0,0.0); macro zero = (0.0,0.0); // -------------------------------------------------------------------------- // module EigenvaluesShiftQRC([N][N](real * real) ?m,a,[N](real * real) e,event !rdy) { [N][N](real * real) q,r; // matrices q,r of QR-decomposition (real * real) trc2,det,dis,rtdis,e1,e2,p; // to determine shift parameter p real dist1,dist2; // distances of eigenvectors for Wilkinson 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 { // determine Wilkinson shift p: the idea is to compute the eigenvalues // of the matrix [[a[N-2][N-2],a[N-2][N-1]],[a[N-1][N-2],a[N-1][N-1]]] // and to select that eigenvalue that is closer to a[N-1][N-1] as shift p // note that this is the only step where complex numbers may be introduced // for real matrices, otherwise the QR-iteration would fail for real matrices // that have complex eigenvalues let(a00 = a[n-2][n-2]) let(a01 = a[n-2][n-1]) let(a10 = a[n-1][n-2]) let(a11 = a[n-1][n-1]) { trc2 = cdiv(cadd(a00,a11),(2.0,0.0)); det = csub(cmul(a00,a11),cmul(a01,a10)); dis = csub(cmul(trc2,trc2),det); // calling croot directly caused inaccuracies, so the following is better rtdis = (dis.1==0.0 ? (dis.0<0.0 ? (0.0,rroot(-dis.0)) : (rroot(dis.0),0.0)) : croot(dis)); e1 = cadd(trc2,rtdis); e2 = csub(trc2,rtdis); dist1 = sqabs(csub(e1,a11)); dist2 = sqabs(csub(e2,a11)); p = (dist1<dist2 ? e1 : e2); } // shift matrix a for(i=0..N-1) next(a[i][i]) = csub(a[i][i],p); pause; // decompose shifted matrix A_k-p*E into factors Q_k and R_k DecompQRC(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 : zero)) let(new_aij_re = Re(p0) + sum(k=0..N-1) Re(cmul(r[i][k],q[k][j]))) let(new_aij_im = Im(p0) + sum(k=0..N-1) Im(cmul(r[i][k],q[k][j]))) next(a[i][j]) = (new_aij_re,new_aij_im); } } pause; } while(sqabs(a[n-1][n-2])>eps*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]) = zero; next(a[n-1][i]) = zero; } pause; } while(n>1); e[n-1] = a[n-1][n-1]; emit(rdy); } drivenby m0 { m[0] = [(1.0,0.0),(2.0,0.0),(3.0,0.0),(4.0,0.0),(5.0,0.0)]; m[1] = [(2.0,0.0),(3.0,0.0),(4.0,0.0),(5.0,0.0),(1.0,0.0)]; m[2] = [(3.0,0.0),(4.0,0.0),(5.0,0.0),(1.0,0.0),(2.0,0.0)]; m[3] = [(4.0,0.0),(5.0,0.0),(1.0,0.0),(2.0,0.0),(3.0,0.0)]; m[4] = [(5.0,0.0),(1.0,0.0),(2.0,0.0),(3.0,0.0),(4.0,0.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,0.0), (-85.0,0.0), (225.0,0.0), (-274.0,0.0), (120.0,0.0)]; m[1] = [ unit, zero, zero, zero, zero]; m[2] = [ zero, unit, zero, zero, zero]; m[3] = [ zero, zero, unit, zero, zero]; m[4] = [ zero, zero, zero, unit, zero]; 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,0.0), (3.0,0.0), (-6.0,0.0), (-2.0,0.0), (4.0,0.0)]; m[1] = [ unit, zero, zero, zero, zero]; m[2] = [ zero, unit, zero, zero, zero]; m[3] = [ zero, zero, unit, zero, zero]; m[4] = [ zero, zero, zero, unit, zero]; 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,0.0), (-14.0,0.0), (-24.0,0.0), (-33.0,0.0), (-18.0,0.0)]; m[1] = [ unit, zero, zero, zero, zero]; m[2] = [ zero, unit, zero, zero, zero]; m[3] = [ zero, zero, unit, zero, zero]; m[4] = [ zero, zero, zero, unit, zero]; 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,0.0), (-2.0,0.0), (4.0,0.0), (3.0,0.0), (-6.0,0.0)]; m[1] = [ unit, zero, zero, zero, zero]; m[2] = [ zero, unit, zero, zero, zero]; m[3] = [ zero, zero, unit, zero, zero]; m[4] = [ zero, zero, zero, unit, zero]; await(rdy); }