Computes the LU factorization on banded matrix A Input arguments: A, square banded matrix band, number of sub and super diagonals which are non-zero Output arguments: L, square matrix of the same dimensions as A containing the lower triangular portion of the LU factorization U, square matrix of the same dimensions as A containing the upper triangular portion of the LU factorization
0001 function [L,U] = LU_banded( A,band ) 0002 % Computes the LU factorization on banded matrix A 0003 % Input arguments: 0004 % A, square banded matrix 0005 % band, number of sub and super diagonals which are non-zero 0006 % Output arguments: 0007 % L, square matrix of the same dimensions as A containing the lower 0008 % triangular portion of the LU factorization 0009 % U, square matrix of the same dimensions as A containing the upper 0010 % triangular portion of the LU factorization 0011 0012 % check user inputs 0013 [m,n]=size(A); % find the size of A 0014 if n~=m 0015 error('input is not a square matrix'); 0016 end 0017 if ~isbanded(A,band,band) 0018 error('matrix is not banded') 0019 end 0020 b=band; 0021 if b>n-1 0022 b=n-1; 0023 end 0024 0025 L=eye(n); % initializing L to the identity matrix 0026 U=A; % initializing U to be A 0027 for k=1:n % loop calculates one column of L and one row of U at a time 0028 % Note that U holds in its lower portion a modified portion of A 0029 for j=k+1:min([k+b,n]) 0030 % if U(k,k) = 0, do nothing, because L is already initialized 0031 % to the identity matrix and thus the k-th column is the k-th 0032 % standard basis vector 0033 if abs(U(k,k)) > 1e-12 % not comparing to zero because of 0034 % possible rounding errors 0035 L(j,k)=U(j,k)/U(k,k); % let the k-th column of L be the k-th 0036 % column of the current U scaled by 0037 % the diagonal element 0038 end 0039 U(j,:)=U(j,:)-L(j,k)*U(k,:); % adjust U by subtracting the 0040 % outer product of of the k-th 0041 % column of L and the k-th row 0042 % of U 0043 end 0044 0045 end