function [L,U,P]=nma_LU(A,threshold) %function [L,U,P]=nma_LU(A,threshold) % %does LU decomposition with permutation matrix for %pivoting reorder. Supports threasold parameter. %Similar to matlab lu function, but with little more control to the %user as to when row exhanges should be made. % %INPUT: % A: an nxn square matrix % threshold: % numerical positive value. row exchanges will be made % only if the abs difference between the largest pivot and % the current pivot is larger than this threshold. % Hence setting threashold to be 0 will cause a row exchange % anytime when there is a larger pivot. This is the DEFAULT % behaviour similar to Matlab lu(). %OUTPUT: % P: nxn permutation matrix such that PA=LU % L: unity lower triangular matrix % U: unity upper triangular % By Nasser M. Abbasi % HW 4, MATH 501, CSUF % % LU decomposition with threshold support. % Copyright (C) 2007 Nasser Abbasi % % This program is free software; you can redistribute it and/or % modify it under the terms of the GNU General Public License % as published by the Free Software Foundation; either version 2 % of the License, or (at your option) any later version. % % This program is distributed in the hope that it will be useful, % but WITHOUT ANY WARRANTY; without even the implied warranty of % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the % GNU General Public License for more details. % % You should have received a copy of the GNU General Public License % along with this program; if not, write to the Free Software % Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. %EXAMPLE RUN % % >> A=rand(4); % >> [L,U,P]=nma_LU(A,0) % if nargin ~=2 error '2 parameters are expected' end if ~isnumeric(A) error 'input matrix A must be numeric' end if ~isnumeric(threshold) error 'threshold must be numeric' end if threshold<0 error 'threshold must be positive' end [nRow,nCol]=size(A); if nRow ~= nCol error 'Matrix must be square' end %****************************** %* Internal functions to flip %* rows for row pivoting %****************************** function flipRows() [c,I]=max(abs(A(n:end,n))); I=I+(n-1); tmp=A(n,:); A(n,:)=A(I,:); A(I,:)=tmp; %make sure we also flip the L matrix rows to keep in sync tmp=L(n,:); L(n,:)=L(I,:); L(I,:)=tmp; %now make the elementary matrix for this move E(n,:)=0; E(n,I)=1; E(I,:)=0; E(I,n)=1; end P=diag(ones(nRow,1)); U=zeros(nRow); L=zeros(nRow); for n=1:nRow-1 currentPivot=A(n,n); E=diag(ones(nRow,1)); maxPivot=max(A(n+1:end,n)); if abs(currentPivot)=threshold flipRows(); end end end P=P*E; %update the perumtation matrix for i=n+1:nRow L(i,n)=A(i,n)/A(n,n); A(i,n)=0; for j=n+1:nRow A(i,j)=A(i,j)-L(i,n)*A(n,j); end end end L=L+diag(ones(nRow,1)); P=P'; U=A; % that is all end