// ************************************************************************** // // // // 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 // // // // ************************************************************************** // // This module computes a QR-decomposition of a given NxN matrix by means of // // Givens rotations. In contrast to module DecompQR that works with reals, the// // version below deals with complex numbers. // // ************************************************************************** // macro N = 5; // -------------------------------------------------------------------------- // // numeric precision // -------------------------------------------------------------------------- // macro eps = 1.0e-8; // -------------------------------------------------------------------------- // // macros on real numbers // -------------------------------------------------------------------------- // macro rroot(x) = exp(x,0.5); macro rsign(x) = (x<0 ? -1 : +1); // -------------------------------------------------------------------------- // // 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.0*x.0 + x.1*x.1); macro cabs(x) = rroot(x.0*x.0 + x.1*x.1); 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 DecompQRC([N][N](real * real) ?m,q,r,event !rdy) { (real * real) a,b,Rq,R,s,c; nat{N} i; // initially, q is the unit matrix and r is m for(i=0..N-1) for(j=0..N-1) { q[i][j] = (i==j ? unit : zero); r[i][j] = m[i][j]; } pause; // next apply Givens rotations to zero out sub-diagonal elements of r for(k=0..N-2) // columns for(jup=0..N-2-k) let(j=N-1-jup) { // rows j=N-1,N-2,…,k+1 if(r[j][k]!=zero) { i = j-1; a = r[i][k]; b = r[j][k]; Rq = cadd(cmul(a,a),cmul(b,b)); R = (Rq.1==0 ? (Rq.0>=0.0 ? (rroot(Rq.0),0.0) : (0.0,rroot(Rq.0))) : croot(Rq)); c = cdiv(a,R); s = cdiv(b,R); // apply Givens rotation to q and r for(l=0..N-1) { next(q[i][l]) = cadd(cmul(c,q[i][l]),cmul(s,q[j][l])); next(q[j][l]) = csub(cmul(c,q[j][l]),cmul(s,q[i][l])); next(r[i][l]) = cadd(cmul(c,r[i][l]),cmul(s,r[j][l])); next(r[j][l]) = (l==k ? zero : csub(cmul(c,r[j][l]),cmul(s,r[i][l]))); } pause; } } // finally, compute the transpose of q for(i=0..N-1) for(j=0..N-1) next(q[i][j]) = q[j][i]; pause; emit(rdy); } drivenby { 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 { 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); }