lu_pp_sl_fl.m

来自「DEMO_COND demonstrates the role of the 」· M 代码 · 共 91 行

function [b, iflag] = lu_pp_sl_fl( A, b, ipivt, m, arith );
% 
% LU_PP_SL   computes the solution of an N by N linear
% system. The LU--decomposition with partial pivoting
% of the system matrix A has to be computed by LU_PP_FL
% before LU_PP_SL_FL is called. 
%
% USAGE:
%
% function [b, iflag] = lu_pp_sl_fl( A, b, ipivt, m, arith )
%
% INPUT:
%    A:     the LU-decomposition of  A  computed by lupiv
%
%    ipivt: pivot information for the LU-decomposition of  A  
%           computed by lupiv
%   
%    b:     the right hand side b
%
%    m:     integer
%           mantissa length
% 
%    arith  string
%           arith = 'c'   chopped arithmetic
%           arith = 'r'   rounded arithmetic
%
%
% OUTPUT:
%    b:     the solution of the linear system if iflag = 0
%
%    iflag: error flag
%           iflag = 0  The system could be solved. The solution
%                      is returned in b.
%           iflag = 1  dimensions of A and b do not match.
%           iflag > 1  zero diagonal element of U
%                      detected in row iflag+1
%
%  Matthias Heinkenschloss
%  Department of Computational and Applied Mathematics
%  Rice University
%  Feb 22, 2001
%
%



iflag = 0;

% get size of A and check dimensions
[nr,nc]    = size(A);
if ( nr ~= nc | nr ~= size(b,1) )
   iflag = 1;
   return
end
n = nr;

% Solve L y = P b  ( b is overwritten by the solution )
 
for k = 1:n-1

%  compute P_k b 
   l = ipivt(k);
   if( k ~= l ) 
       tmp  = b(k);
       b(k) = b(l);
       b(l) = tmp;
   end

%  compute M_k b  
   for i = k+1:n
       tmp  = sflop( A(i,k), b(k), m, arith, '*' );
       b(i) = sflop( b(i), tmp, m, arith, '+' );
   end
end


% Solve U x = y  ( b is overwritten by the solution )
 
for k = n:-1:1
    if ( A(k,k) == 0 )
         iflag = k+1;
         return
    end
    b(k) = sflop( b(k), A(k,k), m, arith, '/' );
    for i = 1:k-1
        tmp  = sflop( A(i,k), b(k), m, arith, '*' );
        b(i) = sflop( b(i), tmp, m, arith, '-' );
    end
end