// ************************************************************************** // // // // 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 computes a LU decomposition of a given NxN matrix as // // described in Cormen, Leiserson, Rivest and Stein's textbook on algorithms. // // Hence, for the given matrix a, the matrices l and u are computed so that // // the equation a = l*u holds, l is a lower left triangular, and u is a upper // // right triangular matrix. // // ************************************************************************** // macro N = 4; module DecompLU([N][N]real ?a,b,l,u,event !rdy) { // store input matrix a in matrix b for(i=0..N-1) for(j=0..N-1) b[i][j] = a[i][j]; // compute LU decomposition of b for(k=0..N-1) { let(piv=b[k][k]) { u[k][k] = piv; for(i=k+1..N-1) { l[i][k] = b[i][k]/piv; u[k][i] = b[k][i]; } } for(i=k+1..N-1) { for(j=k+1..N-1) next(b[i][j]) = b[i][j] - l[i][k] * u[k][j]; } cp: pause; } for(i=0..N-1) l[i][i] = 1.0; emit(rdy); } drivenby { a[0][0] = 2.0; a[0][1] = 3.0; a[0][2] = 1.0; a[0][3] = 5.0; a[1][0] = 6.0; a[1][1] = 13.0; a[1][2] = 5.0; a[1][3] = 19.0; a[2][0] = 2.0; a[2][1] = 19.0; a[2][2] = 10.0; a[2][3] = 23.0; a[3][0] = 4.0; a[3][1] = 10.0; a[3][2] = 11.0; a[3][3] = 31.0; await(rdy); // check the result, i.e. that a = L * U holds, where L and U are stored in b for(i=0..N-1) for(j=0..N-1) assert(a[i][j] == sum(k=0..N-1) (l[i][k] * u[k][j])); }