function x=nma_pentaSolve(A,b) %function x=nma_pentaSolve(A,b) %Solves a penta-diagonal Ax=b system % %INPUT: % A an nxn Matrix % b an vector of length n % %OUTPUT % x vector of length n, the solution for Ax=b % %Algorithm Overview: % The matrix is banded matrix. Due to this, saving in processing % is achived by only processing the non-zero elements along the band. % see full report for more details %by Nasser M. Abbasi. Feb 27,2007 %March 7, 2007 if nargin ~=2 error('2 arguments required'); end if ~isnumeric(A)|~isnumeric(b) error('A and b must be numeric'); end [nRow,nCol]=size(A); if nRow~=nCol error('A must be square'); end b=b(:); [b_nRow,b_nCol]=size(b); if b_nCol>1 error('b must be a vector'); end if b_nRow~=nRow error('b must be same size as A'); end [U,b_new]=penta_elimination(A,b); x=nma_pentaBackSub(U,b_new); end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % The algorithm checks for the band width, and process elements % within the band as using standard G.E. This saves processing % time compared with the G.E. which process the whole matrix. % % Since the matrix is sparse, we save time by processing only the % non-zero band. % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function [A,b]=penta_elimination(A,b) TRUE=1; FALSE=0; [nRow,nCol]=size(A); pivot=1; more_bands=TRUE; while more_bands r=find(abs(A(pivot:end,pivot))>4*eps); r=r+(pivot-1); %find is relative to current band, so adjust to A end_row=r(end); r=find(abs(A(end,:))>4*eps); end_column=r(end); [U,b_new]=nma_gaussian_elimination... (A(pivot:end_row,pivot:end_column),b(pivot:end_row)); A(pivot:end_row,pivot:end_column)=U; b(pivot:end_row)=b_new; if pivot==nRow more_bands=FALSE; else pivot=pivot+1; end end end