This is my naive LU decomposition routine: >> help LU [L,U] = LU(A) LU takes a square matrix A and returns upper and lower triangular matrices such that LU = A. NOTE: it writes over A! This is my naive upper triangular routine: >> help up_solve [x] = up_solve(U,b) takes an upper triangular matrix U and vector b and solves Ux=b This is my naive lower triangular routine: >> help down_solve [x] = down_solve(L,b) takes a lower triangular matrix L and vector b and solves Lx=b Define a random problem >> A = rand(100,100); b = rand(100,1); >> A_o = A; >> flops(0) >> t = cputime; >> [L,U] = LU(A); test the time it took to run in CPU: >> cputime-t ans = 14.0500 test the number of flops: >> flops ans = 666700 >> [y] = down_solve(L,b); >> [x] = up_solve(U,y); My residual is on the order of 10^-12 >> max(abs(A_o*x-b)) ans = 2.0985e-12 Reinitialize A >> A = A_o; Here is matlab's LU decomposition routine. Notice something _really_ annoying about their help page: they refer to the routine as "LU" in the following. But if you were to type [L,U] = LU(A), either you'd get an error, or you'd get my LU.m. (If you'd saved it in the right place.) _All_ matlab routines are in lower-case. >> help lu LU LU factorization. [L,U] = LU(X) stores a upper triangular matrix in U and a "psychologically lower triangular matrix" (i.e. a product of lower triangular and permutation matrices) in L, so that X = L*U. X must be square. [L,U,P] = LU(X) returns lower triangular matrix L, upper triangular matrix U, and permutation matrix P so that P*X = L*U. LU(X), with one output argument, returns the output from LINPACK'S ZGEFA routine. LU(X,THRESH) controls pivoting in sparse matrices, where THRESH is a pivot threshold in [0,1]. Pivoting occurs when the diagonal entry in a column has magnitude less than THRESH times the magnitude of any sub-diagonal entry in that column. THRESH = 0 forces diagonal pivoting. THRESH = 1 is the default. See also LUINC, QR, RREF. >> flops(0) >> t = cputime; >> [L,U] = lu(A); >> cputime-t ans = 0.0400 >> flops ans = 661650 note: Matlab's LU gives back L with pivoting permutations. And rather than untangling all that, I'll just crudley use matlab's matrix inversion routine "inv" to solve the problems. >> y = inv(L)*b; >> x = inv(U)*y; >> max(abs(A_o*x-b)) ans = 7.4774e-14 Note: matlab's LU and inv spank my naive versions. Their residual is 100 times smaller (2.1e-12 vs 7.5e-14). Also, it's got a lower flop count (666700 vs 661650). The marked difference is in the CPU time: matlab's version is over 350 times faster(14.05 vs .040). This is for very good reasons that either I or Will Dickinson can explain. >> diary off