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;   returnendn = nr;% Start LU--decompositionfor 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   endend