// ************************************************************************** // // // // 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 implements the Gram-Schmidt algorithm to compute an // // orthonormal basis of a vector space given any basis of linearly independent// // vectors as columns of the input matrix a. The result is also given in the // // columns of the output matrix b. It is well-known that the algorithm is not // // that good in terms of numeric accuracy. // // ************************************************************************** // macro M = 4; macro N = 3; macro equal(a,b) = (-0.0000001 <= a-b) & (a-b<=0.0000001); macro sqrt(x) = exp(x,0.5); module GramSchmidt([M][N]real ?a,c,event !rdy) { [M][N]real b; // to store the input matrix [N]real ip; // to store the inner products // store input matrix a in matrix b for(i=0..M-1) for(j=0..N-1) b[i][j] = a[i][j]; // start with the first vector as normalized vector of column 0 in b let(length = sqrt(sum(j=0..M-1) (b[j][0] * b[j][0]))) for(i=0..M-1) next(c[i][0]) = b[i][0]/length; pause; for(k=1..N-1) { // compute inner products of b[*][k] with existing orthogonal vectors for(i=0..k-1) ip[i] = sum(j=0..M-1) (b[j][k] * c[j][i]); // generate new orthogonal vector c[*][k] without normalizing for(i=0..M-1) c[i][k] = b[i][k] - sum(j=0..k-1) (ip[j] * c[i][j]); // normalize the computed vector c[*][k] in the next step let(length = sqrt(sum(j=0..M-1) (c[j][k] * c[j][k]))) for(i=0..M-1) next(c[i][k]) = c[i][k]/length; pause; } emit(rdy); } drivenby { // first vector a[0][0] = 2.0; a[1][0] = 1.0; a[2][0] = 0.0; a[3][0] = 3.0; // second vector a[0][1] = 4.0; a[1][1] = 2.0; a[2][1] = 1.0; a[3][1] = -1.0; // third vector a[0][2] = 1.0; a[1][2] = 0.0; a[2][2] = 2.0; a[3][2] = -13.0; await(rdy); // check that column vectors of c are orthogonal to each other for(i=0..N-1) { assert(equal(1.0,sum(k=0..M-1) (c[k][i] * c[k][i]))); for(j=i+1..N-1) assert(equal(0.0,sum(k=0..M-1) (c[k][i] * c[k][j]))); } }