// ************************************************************************** // // // // 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. Alternative methods are based on Householder reductions // // or the Gram-Schmidt procedure to determine orthogonal matrices. In general,// // every matrix M can be written as M = Q*R with a orthogonal matrix Q and a // // upper right triangular matrix R. Since Q is orthogonal, all its columns are// // orthogonal to each other and their lengths are 1. Therefore, the transpose // // of Q is the inverse of Q. // // A Givens rotation on i,j with angle phi is a linear function f(x) = G*x // // that maps a N-vector x to a N-vector f(x) such that x is rotated in the // // (i,j)-plane with angle phi. The rotation can be described by a NxN matrix // // G that is the unit matrix with the following four exceptions: // // * G[i][i] = +cos(phi) // // * G[i][j] = +sin(phi) // // * G[j][i] = -sin(phi) // // * G[j][j] = +cos(phi). // // // // These matrices are also called Jacobi matrices. The vector y := f(x) is // // therefore only changed at positions i and j as follows: // // * y[i] = cos(phi)*x[i] + sin(phi)*x[j] // // * y[j] = cos(phi)*x[i] - sin(phi)*x[j] // // * y[k] = x[k] for k!=i and k!=j // // // // It is not necessary to compute phi and then cos(phi) and sin(phi), instead // // it is sufficient to directly compute cos(phi) and sin(phi). If the goal is // // to zero out an element b of a vector, then we can reduce the consideration // // to two rows: // // // // | c s| |a| |R| // // | | * | | = | | // // |-s c| |b| |0| // // // // Clearly, R = sqrt(a*a+b*b) is the length of the vector (a,b) and we have // // c = cos(phi) = a/R and s = sin(phi) = b/R. // // // // A QR-decomposition of a given matrix m is then obtained by Givens rotations// // R := G_0 * G_1 * ... * G_L * m such that all sub-diagonal elements of m are// // zeroed out by the rotations. The matrix Q is the inverse of (G_0 *…*G_L), // // which can be computed as Q = transpose(G_0 * ... * G_L). // // ************************************************************************** // macro N = 5; macro sqrt(x) = exp(x,0.5); macro rabs(x) = (x<0 ? -x : +x); macro eps = 1.0e-8; module DecompQR([N][N]real ?m,q,r,event !rdy) { real a,b,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 ? 1.0 : 0.0); 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]!=0.0) { i = j-1; a = r[i][k]; b = r[j][k]; R = sqrt(a*a+b*b); c = a/R; s = b/R; // apply Givens rotation to q and r for(l=0..N-1) { next(q[i][l]) = c * q[i][l] + s * q[j][l]; next(q[j][l]) = c * q[j][l] - s * q[i][l]; next(r[i][l]) = c * r[i][l] + s * r[j][l]; next(r[j][l]) = (l==k ? 0.0 : c * r[j][l] - 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,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); // check correctness: for(i=0..N-1) for(j=0..N-1) assert(rabs(m[i][j]-sum(k=0..N-1)(q[i][k]*r[k][j]))<eps); } drivenby { 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); // check correctness: for(i=0..N-1) for(j=0..N-1) assert(rabs(m[i][j]-sum(k=0..N-1)(q[i][k]*r[k][j]))<eps); }