/* return : P; A -> L\U */ def lu(A) { N = size(A)[0]; P = vector(N); for ( I = 0; I < N; I++ ) P[I] = I; for ( J = 0; J < N; J++ ) { for ( I = J; I < N; I++ ) if ( A[I][J] != 0 ) break; if ( I == N ) return 0; if ( I != J ) { printf("Exchange %a and %a rows\n",I,J); for ( K = 0; K < N; K++ ) { T = A[J][K]; A[J][K] = A[I][K]; A[I][K] = T; } T = P[J]; P[J] = P[I]; P[I] = T; } for ( I = J+1; I < N; I++ ) { M = A[I][J]/A[J][J]; for ( K = J; K < N; K++ ) A[I][K] -= M*A[J][K]; A[I][J] = M; } printf("J(row)=%a\n%a\n",J,A); } return P; } AL=newmat(3,3,[[1,0,0],[2,1,0],[3,2,1]])$ AU=newmat(3,3,[[2,-2,1],[0,3,-3],[0,0,1]])$ A=AL*AU; A0=A*1$ printf("P=%a\n",lu(A))$ print(A)$ printf("U is stored in the diagonal and upper part of A.\nL is diag(1)+lower part of A\n")$ printf("Each step is the procedure of making the upper triangular matrix.\n")$ def ex1() { AL=newmat(3,3,[[1,0,0],[2,1,0],[3,2,1]])$ AU=newmat(3,3,[[2,-2,1],[0,0,-3],[0,3,1]])$ A=AL*AU; printf("%a\n",A); A0=A*1$ printf("P=%a\n",lu(A))$ print(A)$ printf("U is stored in the diagonal and upper part of A.\nL is diag(1)+lower part of A\n")$ printf("Each step is the procedure of making the upper triangular matrix.\n")$ } /* A= L\U */ def solve_by_lu(A,P,B) { N = length(B); B1 = vector(N); for ( I = 0; I < N; I++ ) B1[I] = B[P[I]]; for ( J = 0; J < N; J++ ) { for ( K = 0, S = 0; K < J; K++ ) S += A[J][K]*B1[K]; B1[J] -= S; } for ( J = N-1; J >= 0; J-- ) { for ( K = J+1, S = 0; K < N; K++ ) S += A[J][K]*B1[K]; B1[J] = (B1[J]-S)/A[J][J]; } return B1; } end$