// ************************************************************************** // // // // 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 the inverse of a regular NxN matrix a in // // its output matrix b by means of the Gaussian algorithm. // // ************************************************************************** // macro N = 3; module InverseGaussJordan([N][N]real ?a,b,event !rdy) { [N][N]real c; // store input matrix a in matrix c and create 1-matrix in b for(i=0..N-1) for(j=0..N-1) { b[i][j] = (i==j?1.0:0.0); c[i][j] = a[i][j]; } // perform Gaussian-Jordan algorithm for(k=0..N-1) { for(i=0..N-1) { let(lambda = c[i][k]/c[k][k]) { for(j=0..N-1) if(i==k) { next(b[i][j]) = b[i][j]/c[k][k]; next(c[i][j]) = c[i][j]/c[k][k]; } else { next(b[i][j]) = b[i][j] - lambda * b[k][j]; next(c[i][j]) = c[i][j] - lambda * c[k][j]; } } } pause; } emit(rdy); } drivenby { a[0][0] = 1.0; a[0][1] = 2.0; a[0][2] = 0.0; a[1][0] = 2.0; a[1][1] = 3.0; a[1][2] = 0.0; a[2][0] = 3.0; a[2][1] = 4.0; a[2][2] = 1.0; await(rdy); // check the result, i.e. that a * b is the unit matrix for(i=0..N-1) for(j=0..N-1) assert((i==j?1.0:0.0) == sum(k=0..N-1) (a[i][k] * b[k][j])); }