Home > K25104 > LinearSystems > LU_banded.m

LU_banded

PURPOSE ^

Computes the LU factorization on banded matrix A

SYNOPSIS ^

function [L,U] = LU_banded( A,band )

DESCRIPTION ^

 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

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

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

Generated on Mon 18-Jan-2016 10:25:49 by m2html © 2005