tridiagldlsl.m

来自「Interpolation routines in matlab」· M 代码 · 共 47 行

%
% tridiagLDLsl.m
%
% Compute the solution of the symmetric tridiagonal linear system 
%      A x = b
% The LDL^T decomposition without pivoting of the symmetric tridiagonal
% matrix A has to be computed by tridiagLDL.m. 
%
%
% function [b, iflag] = tridiagLDLsl( d, e, b )
%
% input:
%        e:     subdiagonal of L
%
%        d:     diagonal of U
%
% output:
%        b:     solution of the linear system
%
%        iflag: error flag
%               iflag = 0  LU decomposition could be computed
%               iflag = 1  dimension of b does not match dimension of A
%               iflag > 1  zero pivot element detected in row iflag+1
%

function [b, iflag] = tridiagLDLsl( d, e, b );

iflag = 0;

% get size of A and check dimensions
n  = size(d(:),1);
if ( n ~= size(b(:),1) )
   iflag = 1;
   return
end

% Solve  L y = b  (overwite b with result)
for i = 2:n
    b(i) = b(i) - e(i-1) * b(i-1);
end

% Solve  D L^T x = y  (overwite b with result)
b(n) = b(n) / d(n);
for i = n-1:-1:1
    b(i) = b(i)/d(i) - e(i) * b(i+1);
end