lu_sp_fl.m

DEMO_COND demonstrates the role of the condition number of a matrix (with respect to inversion)

function [A, ipivt, iflag] = lu_sp_fl( A, m, arith );
% 
% LU_SP_FL  computes the LU--decomposition with SIMPLE
% pivoting of an N by N  matrix A. 
% All operations are done using m-digit floating point,
% or m-digit chopping arithmetic.

% 
%
% USAGE:
%
% function [A, ipivt, iflag] = lu_sp_fl( A, m, arith )
%
% INPUT:
%    A:     the N by N matrix A
%
%    m:     integer
%           mantissa length
% 
%    arith  string
%           arith = 'c'   chopped arithmetic
%           arith = 'r'   rounded arithmetic
%
%
% OUTPUT:
%    A:     the LU-decomposition of  A
%
%    ipivt: pivot information
%
%    iflag: error flag
%           iflag = 0  Row reductions could be performed,
%                      A is upper triangular.
%           iflag = 1  A is not a square matrix.
%           iflag > 1  zero pivot element 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 )
   iflag = 1;
   return
end
n = nr;

% Start LU--decomposition

for k = 1:n-1

%  Find pivot index, i.e.
%  find the first nonzero entry in column k with row index >= k
   l    = k;
   amax = abs(A(k,k));
   while( amax == 0 & l < n )
       l = l+1;
       amax = abs(A(l,k));
   end
   ipivt(k)   = l;
     

%  Interchange rows if necessary
   if( k ~= l ) 
       tmp      = A(k,k:n);
       A(k,k:n) = A(l,k:n);
       A(l,k:n) = tmp;
   end 

   if( amax > 0 )  % if amax == 0 subcolumn A(k:n,k) is zero
      for i = k+1:n
%          compute multipliers
           A(i,k) = sflop( (-A(i,k)), A(k,k), m, arith, '/' );

%          Perform row elimination
           for j = k+1:n
               tmp    = sflop( A(i,k), A(k,j), m, arith, '*' );
               A(i,j) = sflop( A(i,j), tmp, m, arith, '+' );
           end
      end
   end
end